AD-A223  556 


Naval  Ocean  Systems  Center 

Snn  Dtogo,  CA  92152-5000 


OTIC  FiLE  COPY 


Technical  Document  1772 
June  1990 

Specification  for  a 
Standard  Electromagnetic 
Propagation  Model 


Approved  for  pubftc  r  dtotrt*rtton  It  unlmited. 


90  07  1  9  003 


NAVAL  OCEAN  SYSTEMS  CENTER 

San  Diego,  California  92152-5000 

J.  D.  FONTANA.  CAPT,  USN  R.  M.  HILLYER 

Commmdtr  Technical  Director 


ADMINISTRATIVE  INFORMATION 

The  work  described  in  this  document  was  performed  from  January  1989  to  January  1990 
by  the  Tropospheric  Branch,  Code  543,  Naval  Ocean  Systems  Center,  for  the  Naval  Oceano¬ 
graphic  Office,  Stennis  Space  Center,  Bay  St.  Loui? ,  MS. 


Released  by  Under  authority  of 

H.  V.  Hitney,  Head  J.  H.  Richter,  Head 

Tropospheric  Branch  Ocean  and  Atmospheric 

Sciences  Division 


CONTENTS 


1.0  INTRODUCTION  .  1 

2.0  INPUTS,  OUTPUTS,  AND  LIMITS  .  2 

2.1  Inputs  .  2 

2.2  Outputs  .  3 

2.3  Limits  .  3 

3.0  STANDARD  PROPAGATION  MODEL  .  4 

3.1  Optical  Interference  Region  Models  .  5 

3.1.1  Reflection  Coefficient  Models  .  7 

3.1.2  Antenna  Pattern  Factor  Models  .  9 

3.1.3  Ray  Trace  Models  .  10 

3.1.4  Effective  Earth  Radius  Model  .  12 

3.1.5  Optical  Region  Limits  .  12 

3.2  Diffraction/ Intermediate  Region  Models  .  13 

3.2.1  NOSC  Evaporation  Duct  Model  .  13 

3.2.2  NOSC  Surface-Based  Duct  Model  .  17 

3.2.3  Troposcatter  Region  Model  .  18 

3.3  Standard  Propagation  Model  FORTRAN  Program .  20 

4.0  TEST  CASES  .  23 

5.0  REFERENCES  .  26 

APPENDIX.  Standard  Propagation  Model  Program  Listing  .  A-l 


FIGURES 


3-1.  Two-path  optical  interference  region .  6 

3-2  Ray  trace  variables .  II 

3-3.  Example  of  9.6-GHz  height-gain  curves .  16 

3-4.  Height-gain  curve  for  surface-based  duct  of  aioitrary  height .  18 

3-5.  Geometry  for  troposcatter  loss  calculations .  19 


I 


TABLES 


2-1.  Required  EM  system  inputs  .  2 

2-2.  Required  environmental  inputs  .  2 

4-1.  EM  system  inputs  for  test  cases  .  23 

4-2.  Environmental  test  set  data  .  24 

4-3.  Output  data  for  environment  I  .  24 

4-4.  Output  data  for  environment  2  .  25 

4-5.  Output  data  for  environment  3  .  25 


Hi 


1.0  INTRODUCTION 


standard  electromagnetic  (EM)  propagation  model  has  been  developed  at  the  Naval  Ocean 
Systems  Center.  It  provides  the  user  with  a  method  of  assessing  EM  propagation  from  100  MHz  to 
20  GHz  in  the  marine  environment  for  a  variety  of  atmospheric  conditions.  The  software  implemen¬ 
tation  of  the  model  returns  the  pattern  propagation  factor  in  decibels  when  a  user  supplies  the  model 
with  the  proper  environmental  and  EM  system  inputs.  The  pattern  propagation  factor  is  defined  as 
the  ratio  of  the  actual  electric  field  at  some  point  to  the  field  that  would  exist  at  that  point  under 
free-space  propagation  conditions.  Free-space  propagation  is  the  propagation  of  energy  that  would 
occur  if  an  omnidirectional  point-radiating  source  were  placed  in  outer  space.  The  radiated  energy 
would  travel  outward  in  all  directions,  the  wave  fronts  propagating  away  from  the  source  with  the 
same  velocity  in  all  directions.  Obviously  these  conditions  would  not  be  satisfied  if  the  point  source 
was  placed  in  the  near-earth  environment.  Refraction  by  the  atmosphere  ensure:  that  the  energy  is 
not  propagated  with  the  same  velocity  in  all  directions,  and  the  surface  of  the  earth  can  intercept  and 
reflect  some  portion  of  the  energy.  If  the  functional  form  of  the  pattern  propagation  factor  is  known, 
then  the  propagation  loss  is  known,  since  the  calculation  of  the  free-space  loss  is  simple.  The  pattern 
propagation  factor  is  also  useful  because  it  appears  in  the  radar-range  equation  and  can  be  used 
to  assess  radar  system  performancc/The  model  is  described  in  this  report  in  detail  and  an  ANSI 
Fortran  program  using  the  model  is  provided.  The  program  that  implements  the  model  can  be  either 
incorporated  into  an  application  model  that  requires  EM  propagation  information  or  used  as  a 
stand-alone  program. 

../  I  ,  . 
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2.0  INPUTS,  OUTPUTS,  AND  LIMITS 


2.1  INPUTS 

A  number  of  EM  system  and  environmental  inputs  are  required  to  determine  the  pattern  propa¬ 
gation  factor.  The  necessary  EM  system  parameters  are  given  in  Table  2-1.  The  required  environ¬ 
mental  inputs  are  provided  in  Table  2-2.  The  antenna  beamwidth  and  elevation  angle  parameters  of 
Table  2-1  are  not  required  for  an  omnidirectional  antenna  type. 

Table  2-1.  Required  EM  system  inputs. 


Parameter 

Units 

Valid  Input  Range 

Frequency 

MHz 

100.0  to  20,000  0 

Height  of  Transmitting 

Antenna 

m 

1.0  to  100.0 

Radar  Target/  Receiver 

Antenna  Height 

m 

1.0  to  30000.0 

Transmitting  Antenna 

n/a 

Horizontal,  vertical 

Polarization 

or  circular 

Transmitting  Antenna  Type 

n/a 

Omnidirectional,  sin  (x)/jr 

Antenna  Beamwidth 

deg 

cosecant-squared, 
height-finder,  specific 
system  height-finder 
>  0  to  45.0 

Antenna  Elevation  Angle 

deg 

-10.0  to  10.0 

Range 

km 

1  to  1000.0 

Table  2-2.  Required  environmental  inputs. 


Parameter 

Units 

Valid  Input  Range 

Evaporation  Duct  Height 

m 

0.0  to  40.0 

Surface  Wind  Speed 

kt 

0.0  to  50.0 

Height  Array  —  2  to  30 

Elements 

m 

0.0  to  1 0000.0 

M-unit  Array  —  Each  Element 

Corresponding  to  the  Like- 

Number  Height  Array  Element 

M 

0.0  to  2000.0 
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2.2  OUTPUTS 


The  only  output  is  the  calculated  pattern  propagation  factor  (in  decibels)  for  the  specified  inputs 
of  Tables  2-1  and  2-2.  Sample  program  outputs  for  a  variety  of  environmental  and  EM  system 
inputs  are  presented  in  Section  4.0. 

2.3  LIMITS 

The  standard  propagation  model  described  in  this  document  will  return  a  value  of  the  pattern 
propagation  factor  in  decibels  for  EM  system  operational  parameters  within  the  range  of  validity  of 
the  inputs  of  Table  2-1  and  for  environmental  inputs  within  the  range  of  validity  of  Table  2-2. 
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3.0  STANDARD  PROPAGATION  MODEL 


The  simplest  case  of  electromagnetic  wave  propagation  is  the  transmission  of  a  wave  between 
a  transmitter  and  a  receiver  in  free  space.  Free  space  is  defined  as  a  region  whose  properties  are 
isotropic,  homogeneous,  and  loss-free,  i.e.,  away  from  the  influences  of  the  earth's  atmosphere.  In 
free  space,  the  electromagnetic  wave  front  spreads  uniformly  in  all  directions  from  the  transmitter. 

While  the  total  amount  of  energy  transmitted  does  not  vary,  i.e.,  no  losses  to  absorption,  etc., 
the  energy  is  distributed  over  an  ever-enlarging  surface.  Thus  the  energy  level  along  any  one  ray 
decreases  inversely  with  the  square  of  the  sphere’s  radius.  This  is  called  »he  free-space path  loss.  The 
free-space  path  loss  for  isotropic  antennas,  expressed  in  terms  of  frequency,  is 

Lfs  -  32.44  ♦  20  log (r)  +  20  1og(/)  (1) 


for  r  in  kilometers  and /in  megahertz. 


If  nonisotropic  antenna  radiational  patterns  are  considered  within  the  loss  calculations,  the  los 
is  referred  to  as  propagation  loss  rather  than  path  loss.  The  propagation  loss  can  be  described  with 
the  aid  of  the  pattern  propagation  factor,  which  is  defined  as  the  ratio  of  the  actual  field  strength  at 
a  point  in  space  to  the  field  strength  that  would  exist  at  the  same  range  under  free-space  conditions 
with  the  beam  of  the  transmitter  directed  toward  the  point  in  question.  For  simplicity,  the  term 
propagation  factor  is  used  throughout  this  document  to  refer  to  the  pattern  propagation  factor. 
Thus,  the  effects  of  the  transmitter  antenna  pattern  are  implied  in  all  calculations.  Symbolically  the 
propagation  factor,  F,  is  given  by 


F  - 


UM 

\EJ 


(2) 


where  Ea  is  the  magnitude  of  the  electric  field  under  free-space  conditions  and  E  is  the  magnitude  of 
the  field  to  be  investigated  at  the  same  point. 

7  he  propagation  factor  is  a  desirable  quantity,  since  it  is  an  identifiable  parameter  in  most 
radar-detection-range  equations.  It  contains  all  the  information  necessary  to  account  for  such  effects 
as  sea-surface  reflection,  atmospheric  refraction,  scattering  from  inhomogeneities  in  the  atmosphere 
and  diffraction  from  the  bulge  of  the  earth's  surface.  Thus,  if  the  functional  form  of  F  is  known, 
then  the  propagation  loss  at  any  point  can  be  determined,  since  the  calculation  of  the  free-space  field 
is  quite  simple.  The  propagation  loss  (in  decibels),  including  antenna  patterns,  is  equivalent  to 

L  -  Lfs  -  20  log  (F)  (3) 

Three  regions  require  different  methods  for  obtaining  signal  strength  (or,  equivalently,  propaga¬ 
tion  factor  or  loss)  as  a  function  of  range.  The  first  region  is  called  the  optical  interference,  or  opti¬ 
cal,  region  and  extends  roughly  from  the  transmitter  to  the  radio  horizon.  In  the  optical  region, 
propagation  is  dominated  by  two-path  coherent  interference  between  direct  and  surface-reflected 
waves.  The  other  distinct  region  is  the  diffraction/troposcatter  region,  which  begins  just  beyond  the 
radio  horizon.  A  third  region,  called  the  intermediate  region,  lies  between  the  optical  and  the  diffrac¬ 
tion  region.  The  propagation  factor  in  this  region  is  obtained  by  a  linear  interpolation  between 
F  values  in  the  optical  and  diffraction  regions. 
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The  standard  propagation  model  that  will  be  presented  here  assumes  a  single-layer  atmosphere. 
The  assumption  of  a  single-gradient  atmosphere  is  somewhat  restrictive,  since  nature  does  not 
always  provide  such  a  simple  propagation  medium.  The  basic  assumption  of  the  single-layer  model 
is  that  refraction  can  be  treated  by  assuming  that  the  refractive  bending  of  the  EM  rays  can  be 
accounted  for  by  using  an  effective  earth  radius  that  is  different  (usually  larger)  than  the  true  earth 
radius.  Ray  paths  over  such  an  earth  would  then  appear  to  be  straight  lines  rather  than  curved 
paths.  To  obtain  an  equivalent  single-gradient  atmosphere  from  an  arbitrary  one,  a  ray  trace  must 
be  used.  A  ray  is  traced  from  the  transmitter  height  to  some  arbitrarily  distant  point  through  the 
various  atmospheric  layers.  The  height  at  this  range  is  then  used  to  determine  the  equivalent  single¬ 
gradient  atmosphere  that  would  be  required  to  trace  a  ray  to  this  range  and  height.  This  equivalent 
gradient  is  used  to  define  an  effective  earth  radius  factor.  The  procedure  is  explained  in  more  detail 
in  Section  3.1.4. 

In  the  discussion  of  the  models,  all  heights  are  in  meters,  all  ranges  are  in  kilometers,  and  all 
angles  are  in  radians  unless  specifically  stated  otherwise. 

3.1  OPTICAL  INTERFERENCE  REGION  MODELS 

For  naval  EM  systems  operated  near  the  earth's  surface,  the  electric  field  at  a  receiving  antenna 
or  radar  target  is  the  vector  sum  of  the  field  components  which  arrive  at  that  point  via  the  direct  and 
sea-reflected  paths,  as  shown  in  Fig.  3-1.  The  phase  component  of  the  reflected  ray  will  lag  the  phase 
of  the  direct  path  because  of  the  difference  in  path  lengths.  The  total  phase  lag,  0,  is  given  by 

6  =  6  +  <!>  (4) 


where  6  is  the  path-length  difference  and  4>  is  the  phase  change  caused  by  reflection  from  the 
surface.  Here  the  assumption  is  made  that  the  direct  and  sea-reflected  rays  have  very  nearly  the  same 
spatial  direction,  so  that  the  major  factor  in  their  addition  is  the  phase  difference.  Kerr  (1951)  gives 
the  following  expression  for  Fin  the  absence  of  abnormal  absorption  or  refractive  effects: 

F  =  l/(e,)2  +  [f(t2)  D  R)7+2  /3/?/(e,)/(t2)cos(e)}l/2  (5) 

The / («,)  factors  describe  the  (normalized  to  I)  antenna  pattern,  and  the  angles,  are  shown  in 
Fig.  3-1.  D  is  called  the  divergence  factor  and  takes  into  account  the  spherical  nature  of  the  reflect¬ 
ing  surface  R  is  the  reflection  coefficient  of  the  reflecting  surface  (the  ratio  of  the  magnitudes  of  the 
reflected  and  incident  fields).  F  varies  from  maximum  to  minimum  as  the  total  phase  lag,  0,  changes 
by  7 r  and  can  assume  values  between  0  and  2. 

The  path-length  difference,  6 ,  in  radians,  between  the  direct  and  reflected  rays  is  given  by 

6  =  (4.193////  H;  10"5)/r  (6) 

Here  r  is  the  total  ground  range,  and  H,'  and  Hr'  the  effective  antenna  heights.  H,'  and  H,'  are 
shown  in  Fig.  3-1  and  are  given  by 


///=///-  (1000  rf)/(2  ae) 

(m) 

(7) 

///=///-(  1000  f| )/(2  ae) 

(m) 

(8) 

5 


Figure  3-1.  Two-path  optical  interference  region. 


where  H,  and  Hr  are  the  transmitter  and  receiver/target  heights,  respectively.  ae  is  the  effective 
earth  radius  which  is  defined  as  the  effective  earth  radius  factor,  k,  times  the  mean  earth  radius  of 
6371  km.  r(  and  r2  are  the  reflection  point  ranges,  r,  can  be  determined  by  solving  the  cubic 
equation 


2  / !  -  3  r  r  ]  *  [r?  -  0.002  a,(//,  +  Hr)]rt  +  0.002  at  H,  r  -  0  (9) 

This  equation  is  frequently  solved  by  using  a  Newton-method  iterative  technique,  but  also  has  the 
following  formal  solution  when  Hr~>  H, . 

r,  =  r/2  -  p  cos  [(£  ♦  rrV'3]  (10) 


where 


P  =  {(4/3)[0.001  at{H,  ♦  Hr)  *  (r/2)2]} 1/2 


and 


{  =  cos-'  {[0.002  af{  Hr  +  //,)r]/p3} 


(ID 


(12) 
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The  antenna  pattern  factors,/ (</),  require  angular  information  about  the  angles  a  and  0,  as 
shown  in  Fig.  3-1.  The  magnitude,  R,  and  phase  shift,  $,  require  knowledge  of  the  grazing  angle,  d'- 
These  angles,  in  radians,  are 


a  =  0.001  ( H,  -  H,  )/r  -  r/(2«,) 

(13) 

*h  -  0.001  H,'/ri 

(14) 

y  =  r\ /°e 

(15) 

0  -  -y  -  <i> 

(16) 

in  terms  of  the  variables  shown  in  Fig.  3-1.  The  divergence  factor  can  be  calculated  by  using  the 
equation 

D  =  [1  ♦  (2  r,  r:)/(rfl,  dr)]‘/J  (17) 

Because  Eq.  9  only  applies  for  Hr  >  Hr  these  terminal  heights  are  normally  swapped  if  Hr  <  Hr 
for  the  calculation  of  r,  only.  Equation  13  will  give  correct  values  of  a  for  the  antenna  pattern  calcu¬ 
lations  if  the  true  (unswapped)  values  of  H,  and  Hr  are  used  in  this  calculation.  However,  r |  in  Eq. 

IS  must  be  replaced  with  r2  to  obtain  the  correct  antenna  pattern  factor  for  the  reflected  ray  if  the 
terminal  heights  have  been  swapped. 

3.1.1  Reflection  Coefficient  Models 

The  magnitude  and  phase  shift  of  the  reflected  ray  can  be  calculated  as  a  function  of  the  grazing 
angle,  i/r.  The  magnitude,  R,  and  the  phase  shift,  ♦,  of  the  reflected  ray  for  horizontal  and  vertical 
polarizations,  respectively,  are 


R 


i4> 

v 


V 


Rh--  1 

(18) 

*H  =  »T 

(19) 

n2  sin  (dr)  -  [n2 

,  ,1/2 
-  COS*  (d0J 

(20) 

n2  sin  (dr)  ♦  [n2 

,  ,1/2 
-  cos2  (dr)] 

where  n  is  the  (complex)  index  of  refraction  and  the  subscripts  H  and  V  indicate  the  polarization. 
The  reflection  coefficient  for  circular  polarization,  calculated  in  terms  of  the  horizontal  and  vertical 
coefficients,  is 


RC  =  0.5  [R* 
*c  =  4>h  - 


,  1/2 
+  RH  +  2RyRH  cos  (4»w  -  4v)] 

(21) 

sin-7  sin  (4>w  4>  ^)/(2  RC)J 

(22) 
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The  magnitude  of  the  reflected  ray  is  also  affected  by  the  roughness  of  the  reflecting  surface.  Surface 
roughness  is  included  following  the  models  of  Ament  (1953),  Beard  (1961),  and  Barrick  (1971)  by 
using  the  formulas 

R  =  R0  exp  (-2  {[2  n  h  sin  («/r)]/A)2)  (h  *)/A  <0.110  (23) 

R  =  R0  (0.5018913  -  {0.2090248  -  [(/*  t^)/X]  -  0.558 19}2) '/2 

0.110  <  (/i  0)/X  <  0.260  (24) 

R  =  0.15  R0  (h4i)/\  >  0.260  (25) 

where  R0  is  the  reflection  coefficient  for  a  smooth  surface,  h  is  the  root-mean*squared  (rms)  wave 
height,  and  y  is  the  wavelength.  The  rms  wave  height  is  obtained  as  a  function  of  wind  speed  by 
using  the  Phillips  (1966)  ocean-wave  model 

h  =  0.0051  W\  (26) 

for  wind  speed  ( (Pj)  in  m/s. 

The  square  of  the  index  of  refraction  required  to  make  the  calculation  of  R  and  for  vertical 
and  circular  polarizations  is  given  by 

n}  =  t  -  i  (18,000  a) /f  (27) 

where  «  and  a  are  the  ordinary  dielectric  constant  and  conductivity,  respectively,  of  seawater,  and  / 
is  the  EM  system  frequency  in  megahertz.  The  constants  themselves  are  obtained  as  a  function  of 
frequency  by  using  Blake’s  ( 1970)  equations,  as  follows: 

Case  1:  /  <  1500 

e  =  80  (28) 

a  =  4.3  (29) 

Case  2:  1 500  <  /  <  3000 

e  =  80  -  0.00733(/  -  1500)  (30) 

a  -  4.3  ♦  0.00148(/  -  1500)  (31) 

Case  3:  3000  <  /  <  10,000 

€  =  69  -  0.00243(7  -  3000)  (32) 

o  =  6.52  +  0.001314(7  -  3000)  (33) 

For  frequencies  greater  than  10,000  MHz,  the  10,000  MHz  values  are  used. 
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3.1.2  Antenna  Pattern  Factor  Models 


The  remaining  terms  in  Eq.  5 ,/(«,-),  the  normalized  antenna  pattern  factors,  are  determined  as  a 
function  of  the  antenna  pattern  type,  beamwidth,  and  pointing  angle.  Five  different  antenna  types 
can  be  used:  omnidirectional,  sin  (*)/*,  cosecant-squared,  generic  height-finder,  and  specific  system 
height-finder.  The  specific  system  height-finder  antenna  type  is  not  discussed  here.  Antenna  patterns 
for  these  antennas  can  be  implemented  by  replacing  the  antenna  pattern  functions  with  user-supplied 
data  The  simplest  case  is  that  of  the  omnidirectional  antenna  which,  as  its  name  implies,  has  a  gain 
of  unity  in  all  directions.  That  is ,/(m)  =  •  for  all  angles  n. 

The  second  case  is  the  sin  (x)/x  antenna  type.  The  radiation  pattern  of  this  antenna  is  symmetric 
about  the  elevation  (pointing)  angle  of  the  antenna.  The  pattern  factor  for  this  antenna  is  given  by 
Blake  (1970)  as 

f(n)  -  sin(x)/.r  /(m)  >  0.03.  ~Mm0Jt  -•  M  —  Mmax  (34) 

where 

x  =  c  sin  (n  -  m0)  (35) 

and  n0  and  nmax  are  the  elevation  angle  and  maximum  angle  in  the  main  beam,  respectively.  The 
value  of  c  is  chosen  so  that /(/i)  =  0.7071  when  n  -  y.0  ±  BW/2.  where  Alt' is  the  beamwidth.  This 
normalization  ensures  that  the  antenna  half-power  points  {20  log  [/(/u)]  =  -3  dB)  occur  at  /i  =  n0 
±  BW/2,  which  is  the  usual  definition  of  the  beamwidth  of  the  antenna.  That  is 

c  =  1.39157/sin  (B IF/2)  (36) 

Pattern  factor  calculations  are  limited  to  those  angles  within  the  main  beam  of  the  antenna  down  to 
the  -30  dB  level  (/(ji)  >  0.03].  Angles  greater  than 

Mma*  =  M0  ±  tan'1  [-4/(1  +  A)  ,2  (37) 

where  A  =  rr/c,  are  limited  to  a  pattern  factor  of  0.03.  This  is  equivalent  to  an  antenna  with  its  first 
sidelobes  at  -30  dB,  a  condition  easily  achieved  with  modern  antennas. 

The  generic  height-finder  antenna  is  a  special  case  of  the  sin  (x)/x  antenna.  Height-finder  anten¬ 
nas  typically  sweep  the  beam  upward  in  clevatirn.  This  can  be  simulated  by  substituting  the  direct 
ray  angle,  n-  for  the  elevation  angle,  na.  Then /(m)  =  I  for  all  values,  n,  of  the  direct  ray  set.  As  the 
antenna  beam  is  swept  upward,  the  pattern  factor  for  the  reflected  ray  gradually  tapers  to  the  -30  dB 
level. 

A  fourth  antenna  type  is  the  cosecant-squared  antenna.  This  antenna  pattern  is  not  symmetric 
about  the  elevation  angle.  The  pattern  factor  is  calculated  by  using 


/<M)  =  1 

Mo  ^  M  S  Mo  +  BW 

(38) 

/(M)  =  sin  (BIF)/sin  (m) 

M  >  Mo  +  BW 

(39) 

/(M)  =  [1  -  (Mo  -  nVBW] 

/<M)  >  0.03,  M  <  Mo 

(40) 
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This  antenna  pattern  is  different  from  the  sin  (*)/.*  antenna,  since  the  beamwidth  of  this  antenna 
does  not  coincide  with  the  -3  dB,  or  half-power,  points  of  the  antenna.  The  orientation  of  the 
antenna  given  above  is  the  one  that  would  be  used  for  shipboard  radars.  Cosecant-squared  antennas 
used  on  an  airborne  radar  are  normally  oriented  in  the  reverse  sense  so  that  the  first  two  equations 
abt  v.ould  describe  the  direct  ray  angles  below  the  elevation  angle  n0.  The  third  equation  would 
then  dc.cribe  the  beam  taper  above  the  elevation  angle.  The  antenna  orientation  is  not  optional,  and 
the  antenna  is  always  assumed  to  be  that  of  a  surface-based  system. 

3.1.3  Ray  Trace  Models 

The  standard  propagation  model  obtains  the  required  value  of  the  effective  earth  radius  factor, 
k,  by  means  of  a  ray  trace.  The  model  allows  the  user  to  input  an  A/-unit-versus-height  profile, 
which  is  used  in  performing  the  ray  trace.  The  ray  trace  equations  are  based  on  small-angle  approx¬ 
imations  to  Snell's  law  and  on  the  assumption  of  a  linear  variation  of  modified  refractivity,  M,  with 
height  up  to  30  vertical  segments.  The  trace  of  an  individual  ray  begins  with  an  elevation  angle  spec¬ 
ified  at  some  initial  height  and  range  and  consists  of  a  series  of  calculations  to  determine  a  series 
of  height  and  range  points  along  the  ray  trajectory.  The  A/-unit  profile  is  constructed  so  that  the 
Af-unit  value  at  the  surface  and  a  zero-meter  height  are  the  first  elements  in  the  profile  arrays.  The 
remainder  of  the  A/-unit  profile  has  a  height  array  of  ascending  heights,  in  meters,  and  an  Af-unit 
array  with  the  corresponding  A/-unit  value  in  a  like-numbered  array.  A  third  array  can  be  con¬ 
structed  from  these  two  arrays  which  contains  the  gradient  between  adjacent  layers.  The  general 
definition  for  this  array  is 


dMdh ,  =  I0*3(  A/, , ,  -  A/,)/(//,m  -  Ht)  (41) 

where  Af;  denotes  the  Af-unit  array  and  Hj  the  height  array  elements,  respectively.  Negative  values 
of  dMdht  indicate  trapping  layers.  A  standard  atmosphere  (4/3  earth)  gradient  is  usually  defined  for 
the  gradient  above  the  highest  height  array  element,  that  is.  dMdhk  -  0.0001 18,  wheie  k  is  the  index 
of  the  last  element  in  the  H  array.  dMdh  values  of  zero  are  not  allowed,  which  is  equivalent  to  not 
allowing  the  Af-unit  values  of  adjacent  height  values  to  be  equal. 

A  critical  launch  angle  can  be  determined  for  transmitter  heights  within  ducts.  This  critical  angle 
is  defined  as  the  minimum  positive  launch  angle  not  trapped  in  the  duct.  The  positive  critical  angle  is 
given  by 


ac  =  lO-3[2(A/„,  -  Mmin)]'n  ♦  10  5  (42) 

while  the  minimum  negative  critical  angle  is  equal  to  -ac.  Here  M  H,  is  the  Af-unit  value  at  the 
transmitter  and  M min  is  the  mimimum  Af-unit  value  at  some  height  greater  than  Hr  If  »Y(  is  in  a 
duct,  then  the  duct  is  treated  as  a  surface-based  duct  even  if  it  does  not  extend  to  the  surface.  The 
height  w  here  M min  occurs  is  the  height  of  the  surface-based  duct.  H,  must  be  in  the  the  duct  for 
Eq  42  to  be  valid,  though  if  H ,  is  above  the  duct,  -ac  would  define  the  launch  angle  for  a  ray 
tangent  to  the  top  of  the  duct  at  some  range.  Rays  launched  with  angles  ac  >  a  >  -ac  will  be 
trapped  within  the  duct. 

The  general  ray  trace  equations  using  the  H,  M,  and  dMdh  arrays  can  be  divided  into  three 
categories:  rays  with  the  terminal  range  known,  rays  with  the  terminal  height  known,  and  rays  with 
the  terminal  elevation  angle  known.  Figure  3-2  illustrates  a  ray  with  a  positive  launch  Angle,  but  the 
equations  also  apply  to  negative  launch  angles  when  proper  attention  is  paid  to  the  layer  indices  and 
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-  - - Hi 

Figure  3*2.  Ray  trace  variable!. 


sign  of  the  launch  angle.  The  equations  given  apply  only  to  range  and  height  values  within  individual 
layers.  All  heights  are  in  meters  and  ranges  in  kilometers. 

Case  I:  h'  known.  oi*0 

o'  =  [oJ  ♦  0.002  dMdh,(h’  -  h)]U7  (43) 

f  -  r  *  (o'  -  a)/dMdh,  (44) 

Case  2:  /  known,  a  **  0 

o'  =  o  ♦  dMdh,  (r'  -  r)  (45) 

h’  -  h  ♦  (o'2  -  o2)/(0.002  dMdh,)  (46) 

Case  3:  o'  known 

r'  =  r  +  (o'  -  o )/dMdh,  (47) 

h'  -  h  *  (o'2  -  o2)/(0,002  dMdh ,)  (48) 

If  the  radicand  of  Eq.  43  is  negative,  there  is  no  solution  for  the  given  height,  h\  since  the  ray  has 
reached  a  maximum  (or,  in  the  cate  of  a  Jowngoing  ray,  a  minimum)  height  less  (or  greater)  than  h'. 
In  this  cate,  the  range  and  height  of  the  ray  maximum  (minimum)  are  given  by 

r'  ~  r  -  a/dMdh,  (49) 


h'  -  h  -  o2/(0.002  dMdh,) 


(30) 


while  a'  =  0  at  this  range  and  height.  One  unique  case  not  covered  by  the  above  equations  is  the 
special  case  a  =  0.  In  this  case,  if  dMdht  >  0,  the  ray  will  become  an  upgoing  ray;  if  dMdht  <  0,  the 
ray  will  become  a  downgoing  ray.  Equations  43  through  SO  can  be  iteratively  used  to  trace  ray  paths 
through  the  user-specified  stratified  atmosphere. 

3.1.4  Effective  Earth  Radius  Model 

The  ray  trace  equations  are  used  to  trace  a  ray  through  the  specified  atmosphere  to  deter¬ 
mine  the  value  of  the  effective  earth  radius  factor,  k,  since  k  is  only  defined  for  a  single-layer  atmo¬ 
sphere.  A  ray  with  a  launch  angle  a  =  0,  or  a  =  ac  if  the  transmitter  is  in  a  duct,  is  traced  to  370  km 
(200  nmi).  The  ray  height  at  that  range  is  compared  to  the  height  a  ray  would  reach  in  a  standard 
atmosphere  to  calculate  a  value  of  k.  The  single-layer  equivalent  A/-unit  gradient  is  given  by 

Q  =  -(2a)/370  ♦  (0,002  Atf)/(370)2  (51) 

where  AH  is  the  terminal  height  of  the  ray  at  370  km  minus  the  launch  height,  Hr  and  a  is  the  ray 
launch  angle.  Then  k  is  defined  as 

k  -  \ HQ  a)  1  <  k  <  5  (52) 

where  a  is  the  mean  earth  radius  of  6371  km.  In  a  standard  atmosphere,  k  -  4/3. 

3.1.5  Optical  Refion  Limits 

The  expression  for  F  in  Eq.  5  is  valid  for  all  values  of  6  so  that  the  path-length  difference,  6,  is 
greater  than  or  equal  to  one-quarter  wavelength  (rr/2  radians),  or  the  grazing  angle  is  equal  to  a 
limit  given  by  Reed  and  Russell  (1966)  at  which  the  spherical  earth  divergence  factor  becomes 
invalid.  This  limit  is  given  by  the  expression 


*Um  =  ,an'‘  f0  001  M/(2  rr  a,)]'  3  =*  (0.0!957)/(*/)'/3  (53) 

where  i k  is  the  grazing  angle, is  the  wavelength  in  meters,  at  is  the  effective  earth  radius,  and /is 
the  frequency  in  megahertz. 

The  optical  region  maximum  range  may  be  reduced  from  that  calculated  for  the  applicable  opti¬ 
cal  region  limit  above  if  the  surface-based  duct  height  is  not  zero.  If  the  transmitter  is  in  a  surface- 
based  duct,  then  the  end  of  the  optical  region  must  correspond  to  a  direct-ray  elevation  angle  greater 
than  or  equal  to  the  critical  angle,  ac.  If  the  value  of  the  elevation  angle,  a,  corresponding  to  the 
quarter-wavelength  or  grazing  angle  limit  value  of  9,  above,  is  not  greater  than  ac,  then  the  optical 
region  limit  is  taken  to  be  the  first  optical  region  peak  (9  =  2 it,  4n,  etc.)  that  is  associated  with  a 
value  of  a  greater  than  af,  This  range  is  determined  by  finding  the  value  of  9  as  a  function  of  ae. 
9(ar)  can  be  obtained  by  using  the  ray  trace  equations  to  find  the  range  corresponding  to  ac,  that  is 

r(af)  =  (K2  +  0.002 (Hr  -  H,)/at]il2  -  ac)af  (54) 

where  dAfdh,  =  l/o,  for  the  single-gradient  atmosphere.  Once  the  range  is  known,  the  value  of  9  is 
obtained  by  determining,  in  order,  the  reflection  point  range,  r,,  the  effective  terminal  heights,  the 
path-length  difference,  6,  the  grazing  angle,  0,  and  the  magnitude  of  phase  change  caused  by  reflec¬ 
tion,  <t>  If  9(ar)  is  less  than,  or  equal  to,  2rr  (the  first  optical  region  peak)  and  greater  than  the 
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quarter-wavelength  or  grazing  angle  limit,  then  6 (ac)  is  the  optical  region  limit.  If  8(arf)  is  greater 
than  2rr,  then  the  next  optical  region  peak  with  8  =  n2rr  >  8faf),  n  =  2,  3,  4,  etc.,  is  taken  to  be 
the  end  of  the  optical  region.  (This  range  is  determined  by  decreasing  the  range  from  that  of  r(ac) 
and  determining  a  new  value  of  8  for  each  range  until  the  desired  range  is  found.)  If  both  terminal 
heights  are  inside  the  surface-based  duct,  then  the  end  of  the  optical  region  is  defined  by  the  quarter- 
wavelength  or  grazing  angle  limit  as  before. 

3.2  DIFFRACTION/INTERMEDIATE  REGION  MODELS 

Beyond  the  horizon,  the  chief  contributions  to  the  electric  field  are  from  diffraction  and,  at 
somewhat  greater  ranges,  tropospheric  scatter.  The  diffraction  field  can  be  represented  as  a  sum  over 
the  possible  number  of  modes  which  are  the  solution  to  the  fundamental  equation  of  mode  theory. 
For  a  standa'd  atmosphere,  the  series  describing  the  field  converges  rapidly  and  only  a  single  mode 
is  necessary  .o  adequately  determine  the  field.  A  single  mode  may  also  describe  the  field  in  the  pres¬ 
ence  of  evaporation  ducts  or  surface-based  ducts  caused  by  elevated  layers,  especially  the  former. 
However,  close  to  the  horizon,  the  series  solution  converges  rather  slowly.  This  is  the  “intermediate 
region,”  and  a  method  of  "bold  interpolation"  originally  described  by  Kerr  (1951)  is  used  to  estimate 
the  field  in  this  region.  This  method  involves  a  linear  interpolation  on  the  logarithm  of  the  pattern 
propagation  factor  from  the  last  valid  range  in  the  optical  region  to  the  first  range  in  the  diffraction 
region. 

The  minimum  range  at  which  the  diffraction  field  solutions  are  applicable  and  the  intermediate 
region  ends  is  given  by  Reed  and  Russell  (1966)  as 


rd  =  Thor  +  230.2  (k2/f)],i  (km) 


(55) 


where  the  horizon  range  is  given  by 


rhor  =  3.572  [(A  H,)'*  ♦  (k  Hr)'2]  (km)  (56) 

for  Ht  and  Hr  in  meters  A  minimum  effective  earth  radius  of  k  -  1.33  is  assumed  for  the  calculation 
of  rd. 

The  diffraction/intermediate  region  models  are  used  to  determine  propagation  loss  as  a  function 
of  height  and  range  for  ranges  and  heights  below  the  lower  angular  limit  of  the  optical  interference 
region.  There  are  three  models  used  to  calculate  loss  in  this  region.  If  the  surface-based  duct  height  is 
zero,  then  the  loss  is  calculated  by  using  a  model  derived  from  the  NOSC  waveguide  program.  If  a 
surface-based  duct  is  present,  an  empirical  model  is  used  to  calculate  loss.  At  somewhat  greater 
ranges,  troposcatter  loss  is  calculated  by  using  a  model  taken  from  Yeh  (I960).  The  model  has  been 
modified  by  the  addition  of  a  “frequency  gain"  factor  from  Rice,  et  al.  ( 1965)  that  gives  better  values 
for  low-altitude  paths.  The  troposcatter  loss  is  calculated  for  all  ranges  beyond  rd  and  added  to  the 
surface-based  duct  or  evaporation  duct  loss  until  the  troposcatter  loss  is  18  dB  less  than  the  applica¬ 
ble  loss.  Beyond  that  point,  only  the  troposcatter  loss  is  calculated. 

3.2.1  NOSC  Evaporation  Duct  Model 

The  evaporation  duct  loss  (in  decibels)  may  be  written  as 

L  -  51.1  +  r  -  Fzl  -  F„  +  10  log  (p)  +  ap  -  Ld  (57) 
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The  loss  term,  Ld,  is  determined  by  using 


Ld  =  20  log  LT(m)1  (58) 

where  the  antenna  pattern  factor,/(p),  gives  a  measure  of  how  much  energy  is  directed  toward  the 
horizon  and  p  represents  the  lowest  direct  ray  angle  in  the  optical  region.  T  is  the  excitation  factor, 
Fu  and  Ftr  the  height-gain  functions  for  the  EM  system  transmitter  and  radar  target/receiver, 
respectively,  p  the  (scaled)  range,  and  a  the  attenuation  rate.  The  specific  values  of  these  quantities 
arc  obtained  as  functions  of  the  duct  height.  The  functions  which  produce  these  values  are  the  result 
of  curve-fitting  the  various  quantities  to  waveguide  program  solutions.  F  is  obtained  by  substituting 
Eq.  57  into  Eq.  3. 

The  waveguide  solutions  used  to  develop  the  evaporation  duct  model  were  made  at  a  single 
frequency,  9.6  GHz.  The  evaporation  duct  solutions  for  other  frequencies  share  a  family  resem¬ 
blance:  the  height  of  the  duct  which  produces  a  particular  propagation  characteristic  varies  inversely 
with  the  frequency.  This  fact  allows  the  solutions  at  9.6  GHz  to  be  scaled  to  other  frequencies.  All 
actual  ranges  and  heights  are  multiplied  by  the  scale  factors 


Rn  =  4.705  I0-2  fn  (59) 

and 

ZN  -  2.214  I0*3  fW  (60) 

respectively,  to  scale  the  solutions  at  other  frequencies  to  the  9.6-GHz  values.  The  coefficients  ensure 
that  Ry  =  ZN  -  I  when  the  frequency  is  set  equal  to  9600  MHz.  When  these  scale  factors  arc  used, 
the  actual  evaporation  duct,  receiver,  and  transmitter  heights  are  scaled  to  the  9.6-GHz  equivalents, 
and  the  range  is  similarly  changed  to  conform  to  the  9.6-GHz  requirements.  For  example,  the  scaled 
duct  height,  A,  is  equal  to  the  actual  evaporation  duct  height,  i,  times  ZN.  Similarly,  if  r  is  the 
actual  range  and  H,  the  actual  EM  system  transmitter  height,  then  the  scaled  range,  p,  is  RN  times  r 
and  the  scaled  transmitter  height,  zr  is  Z  v  times  H,. 

The  height  gains  expressed  as  a  function  of  scaled  duct  height  are  of  two  different  forms, 
depending  on  whether  or  not  the  duct  height  is  sufficient  to  support  a  well-trapped  mode.  The 
height-gain  function  (in  decibels)  for  scaled  duct  heights  less  than  10.25  meters  may  be  written  as 

F(z)  =  Cl  ♦  C3  zc*  ♦  C5  z  >  1.0  (61) 

where  z  is  the  scaled  height  of  either  the  EM  system  transmitter  or  the  radar  target/receiver.  The  C7 
are  constants  that  are  a  function  of  the  scaled  due:  height.  For  well-trapped  modes  (i.e.,  scaled  duct 
heights  between  10.25  and  23.3  meters),  two  functions  are  necessary  to  obtain  the  height-gains  in 
decibels; 


f\z)  =  Cl  ln[sin  (C2  zC3)]  ♦  C4 

1.0  <  r  <  Zmox 

(62) 

F[z)  =  C5  zC6  ♦  C7 

2  '>  max 

(63) 
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As  before,  the  coefficients,  Ci,  are  determined  from  the  scaled  duct  height,  and  z  is  the  scaled  height 
of  the  EM  system  transmitter  or  radar  target/ receiver.  Z is  calculated  by  using  the  formula 


Zmax  *  4e-03l<*-|00>  +  6  (64) 

where  A  is  the  scaled  duct  height.  The  coefficients  for  scaled  duct  heights  less  than  I0.2S  meters  are 
calculated  by  using  the  following  formulas: 

Cl  =  (-2.2  e*° 244A  ♦  \1)A.12~C2  (65) 

C2  =  [4.062361  I04  -  (A  ♦  4.496l)2]'/2  -  201  .0128  (6t>) 

Ci  -  (-33.9  e-° 5I7A  -  3)4.72'C4  (67) 

C4  =  (1.43012  104  -  (A  ♦  5.32545)2]172  -  119.569  (68) 

C5  =  41  e*04,A  +  61  (69) 


The  coefficients  for  scaled  duct  heights  between  10.25  and  23.3  meters  are  calculated  by  applying  the 
following  formulas: 

Cl  -  -0.1189  A  +  5.5495  (70) 

C2  -  {1.3291  sin  [0.2 1 8( A  -  10) 77]  +  0.2171  In  (A))4.72'CJ  (71) 

C3  =  3/2  (72) 

C4  =  87  -  [313.29  -  (A  -  25.3)2]'/2  (73) 

C5  =  Fmaxl{Zmax-“)  (74) 

C6  =  (Zm„/4.72)(S/Fm,J  (75) 

C7  -  49.4e-01699<‘i',0>  +  30  (76) 

where 

5  =  4.72  Cl  C2  Ci  (Zmax)n  /tan  [C2  (Zm„)cn  (77) 

and 


Fmax  =  CI  In  {sin  [C2  (Zmax)a]}  *  C4  -  C7  (78) 

which  are  necessary  to  make  the  two  functions,  F{z),  and  their  slopes  continuous  about  Zmex.  Using 
these  coefficients  in  the  equations  will  produce  height-gain  curves  which  increase  with  height  for 
scaled  duct  heights  below  10.25  meters.  The  well-trapped  modes  have  an  initial  increase  with  height 
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for  a  limited  range  of  z  near  the  surface,  peak,  and  then  decrease  with  height  to  some  value,  there¬ 
after  displaying  very  little  variation  with  height.  The  minimum  scaled  height  used  for  calculating  the 
height  gains  is  1.0  meter,  and  heights  below  this  are  set  equal  to  this  value. 

Scaled  duct  heights  greater  than  23.3  meters  have  more  than  one  mode  which  can  propagate  in 
the  guide.  The  effect  of  the  multiple  modes  is  to  add  constructively  at  some  range/height  combi¬ 
nations  and  destructively  at  others,  a  condition  similar  to  the  optical  region  interference.  Since  this 
variation  is  not  predictable  without  using  a  waveguide  program,  the  scaled  duct  heights  greater  than 
23.3  meters  are  treated  as  23.3  meter  ducts.  An  example  of  height-gain  curves  for  evaporation  ducts 
is  shown  in  Fig.  3-3. 


Two  factors  from  Eq.  57  remain  to  be  specified,  T  and  a.  The  excitation  factor,  which  is  a 
measure  of  the  relative  strength  of  the  mode,  may  be  obtained,  in  decibels,  by  using 


216.7  ♦  1.5526  A 

A 

< 

3.8 

(79) 

222.6  -  1 . 1 77 1  (A  -  3.8) 

A 

> 

3.8 

(80) 

The  attenuation  rate  in  d B/  km  is 

a  =  {92.516  -  [8608.7593  -  {A  -  20.2663)2] '/2J  (81) 

for  values  of  a  S  0.0009,  which  is  the  lowest  attenuation  rate  used.  It  is  convenient  to  replace  the 
attenuation  rate  term  in  Eq.  57,  op,  with  fir,  where  r  is  the  actual  range  and 

fi  =  a  Rn  (82) 
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The  attenuation  rates  for  the  higher  duct  heights  may  be  several  orders  of  magnitude  smaller  than 
the  standard  diffraction  (zero-meter  duct  height)  rate. 

3.2.2  NOSC  Surface-Based  Duct  Model 

The  NOSC  model  for  a  surface-based  duct  from  elevated  layers  is  not  as  complex  as  the  evapo¬ 
ration  duct  model.  It  is  based  on  an  empirical  fit  of  experimental  data.  The  loss  (in  decibels)  may  be 
written 


L  -  C  -  Ftr  +  20  log  (r)  -  Ld  (83) 

where  Fzr  is  the  height-gain  function  for  the  receiver/target  height  and  Ld  is  defined  in  Eq.  58. 

He.c  C  (in  decibels)  is  given  by  C  -  32.44  +  20  log  (/).  The  attenuation  rate  term  is  not  used  in  this 
model,  and  no  range  or  height  scale  factors  are  used  either.  Similarly,  the  only  height-gain  term  used 
in  Eq.  83  is  the  height  gain  of  the  radar  target/receiver  height.  As  the  “guide”  dimensions  are  nor¬ 
mally  on  the  order  of  hundreds  of  meters,  these  ducts  affect  frequencies  of  100  MHz  and  below, 
unlike  the  evaporation  duct,  which  only  affects  frequencies  greater  than  I  GHz.  This  model  has  the 
disadvantage  of  being  anisotropic  with  choice  of  terminal  heights.  The  height-gain  term  used  in  the 
standard  propagation  model  is  always  calculated  for  the  terminal  height,  specified  as  the  radar 
target/receiver  height.  F  is  obtained  by  substituting  Eq.  83  into  Fq.  3. 

The  height-gain  function  for  the  surface-based  duct  model  is  calculated  as  a  function  of  fre¬ 
quency  and  duct  height  for  any  arbitrary  radar  target/receiver  height,  z.  as  follows: 


Case  1:  IOO</<  150 

Fir  =  -60 (z/D  -  0.5)2 

z/D  <  0.8 

(84) 

Ftr  =  1.14(z/D)<>26  -  10 

z/D  >  0.8 

(85) 

Case  2:  150  </<  350 

F„  =  10  -  200 (z/D  -  0.5)4 

z/D  <  1.0 

(86) 

Fzr  =  7.5  (z/D)-133  -  10 

z/D  >  1.0 

(87) 

Case  3:  /  >  350 

F„  -  10  -  200 (z/D  -0.5)4 

z/D  <  1.0 

(88) 

F„  =  12.5(z/D)-«  -  15 

z/D  >  1.0 

(89) 

Here  D  is  the  duct  height.  An  example  of  the  height-gain  curves  produced  by  these  formulas  is  given 
in  Fig.  3-4.  The  shapes  of  the  height-gain  curves  are  characteristic  of  well-trapped  modes,  as  should 
be  expected  from  a  surface-based  duct. 
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GAIN  (dB) 

Figure  3-4.  Height-gain  curve  for  surface-based  duct  of  arbitrary  height. 

3.2.3  Troposcatter  Region  Model 

At  ranges  sufficiently  greater  than  the  horizon,  scattering  from  irregularities  in  the  troposphere 
begins  to  dominate  the  electric  field.  Yeh  ( I960)  gives  the  troposcatter  loss  in  decibels  as 

L  -  114.9  +  0.08984(r  -  rhor)/k  +  20  log  (r) 

♦  30log(/)  -  0.2  AT,  -  Ld  +  Hc  (90) 

Here  r  is  the  range,  rhor  is  the  horizon  range,  N}  is  the  surface  refractivity  value,  and  H0  is  the 
irequency-gain  function  from  Rice  et  al.  (1965).  Ld  is  defined  in  Eq.  58.  Fis  obtained  by  substituting 
Eq.  90  into  Eq.  3. 

The  frequency-gain  function,  H0>  is  primarily  of  importance  for  low  antenna  heights,  especially 
if  the  system  frequency  is  very  low.  The  procedure  for  obtaining  H0  requires  a  calculation  of  the 
effective  scattering  height,  h0,  which  is  equal  to 

A,  *  (ir»)/(l  ♦  »)»  (km)  (91) 

where  r  is  the  ground  range,  0  is  the  scattering  angle,  as  shown  in  Fig.  3-5,  and  s  is  defined  by 

s  =  U  X  10.0  >  J  >  0.10  (92) 
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Figure  3-5.  Geometry  for  troposcatter  loss  calculations. 


The  angles  from  these  equations  are  given  by 


e  = 

e0  ♦  e,  ♦ 

e2 

(93) 

e0  =  r/ae 

(94) 

0|  =  rjae 

(95) 

e2  =  r2/ae 

(96) 

{ =  e/2 

*  e,  ♦  (H, 

-  H,)/r 

(97) 

x  =  e/2 

+  e2  ♦  {H, 

-  H,)/r 

(98) 

in  terms  of  the  effective  earth  radius,  ae,  the  tangent  ray  ranges,  r,  and  r2,  the  terminal  heights,  H, 
and  Hr,  and  the  total  range,  r,  as  shown  in  Fig.  3-5.  The  tangent  ranges,  r t  and  r2,  are  equal  to 


r,  =  (0.002  at  H,)'* 

(km) 

(99) 

1/2 

r2  =  (0.002  a,  Hr) 

(km) 

(100) 

The  frequency-gain  function  is  then  defined  as 

H0  =  +  H0  >  0.0  (dB)  (101) 
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where 


//,  =  [//„(/?,)  +  H0(^2)]/2  (102) 

If  AH0  i:  greater  than  H] ,  then  H0  is  equal  to  twice  the  value  of  Ht .  The  function  is  calculated 
by  using 


*o(*,)  =  Cl  <*l  ♦  Cz)-4/3 

(103) 

H0(R2)  =  c,  ( R2  +  c2)-«'3 

(104) 

where  R|  and  R2  are  functions  of  the  terminal  heights  and  EM  system  frequency,/,  in  megahertz. 
These  variables  are  calculated  as  follows: 

R |  =  0.0419  fH,  6 

R2  =  0.0419  /  Hr  e 

And  the  terms  c,  and  c2  are  defined  as 

=  16.3  +  13.3  vt 
c2  -  9.40  +  0.16 

The  factor  rj,  must  be  calculated  as  a  function  of  h0  in  the  following  manner: 

tj,  =  0.5696  /i0[l  ♦  (0.031  -  0.00232  Ns  «■  5.67  N]  10*«)exp  (-3.8  h60  10'6)] 

5.0  >  r\s>  0.01  (109) 

The  remaining  term,  A H0,  is  calculated  by  using 

=  6  [0.60  -  log  (»;,)]  log(s)  log  (9)  (dB)  (110) 


(105) 

(106) 

(107) 

(108) 


where  q  i»  given  by 

q  -  R2/(s  /?,)  10.0  >  q  >  0.10  (III) 

The  correction  term  A//0  is  zero  for  =  4.0,  s  =  1.0,  or  q  -  1.0  and  has  a  maximum  value  of 
3.6  dB  for  highly  asymmetrical  paths  when  -  1.0. 

3.3  STANDARD  PROPAGATION  MODEL  FORTRAN  PROGRAM 

The  standard  propagation  model  is  implemented  in  a  program  called  FFACTR.  FFACTR  is 
written  in  ANSI  Fortran  77  with  the  allowable  M1L-STD-I753  extensions.  FFACTR  will  return  a 
single  value  for  the  propagation  factor  in  decibels  for  the  specified  EM  system  and  environmental 
parameters  of  Table  2-1  and  Table  2-2.  To  use  FFACTR,  the  operator  must  compile  and  link  the 
routines  that  comprise  the  program.  A  complete  list  of  all  subroutines  is  included  in  the  appendix. 
The  subroutines  are  listed  in  alphabetical  order  following  lists  of  the  MAIN  and  FFACTR  routines 
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and  the  common  block  “include"  files.  No  EM  system  or  environmental  libraries  are  supplied  with 
FFACTR,  though  a  limited  number  of  environmental  and  EM  system  data  files  are  supplied  for  test 
purposes. 

Because  the  standard  propagation  model  is  designed  to  return  a  single  value  for  a  collection  of 
input  data.  FFACTR  behaves  more  like  a  subroutine  than  a  stand-alone  program.  A  demonstration 
program,  MAIN,  which  acts  as  a  driver  for  the  FFACTR  program,  is  included  to  demonstrate  one 
possible  use  of  the  program.  The  driver  simulates  an  inbound  radar  target  by  supplying  FFACTR 
with  constant  environmental  and,  except  for  the  range,  EM  system  parameters.  The  range  decreases 
from  IOC  km  to  10  km  in  10-km  increments.  In  addition  to  MAIN,  several  subroutines  are  provided 
which  allow  the  operator  to  enter  environmental  and  EM  system  data  from  the  keyboard  or  files. 
These  subroutines,  SYSFIL,  ENVFIL,  ENVINP,  and  SYSINP,  are  not  intended  as  part  of  the 
FFACTR  program,  but  only  for  use  in  verifying  the  correct  operation  of  the  FFACTR  program. 

The  operational  sequence  of  FFACTR  is  detailed  in  the  following  paragraphs. 

First,  the  environmental  and  EM  system  data  are  entered  and  the  FFACTR  subroutine  is  called. 
FFACTR  then  calls  the  MPROF  subroutine  to  process  the  AZ-unit,  (A/,),  and  height,  (//,),  array 
data  and  to  construct  the  dMdht  array.  The  profile  is  inspected  for  surface-based  ducts  by  the 
DUCTS  subroutine  and,  if  one  exists  that  contains  the  EM  system  transmitter  height,  the  critical 
angle,  a(,  is  calculated.  Subroutines  PUSH  and  1NSRT  are  used  by  MPROF  to  insert  an  A/-unit 
array  value  at  the  transmitter  height.  Upon  return  from  the  MPROF  subroutine,  the  GETK  subrou¬ 
tine  is  called  to  perform  the  raytrace  that  is  used  to  obtain  the  effective  earth  radius  factor.  Next, 
several  subroutines  are  called  to  initialize  various  constants  required  by  the  program.  ANTPAR 
converts  the  EM  system  antenna  parameters  from  degrees  to  radians  and  establishes  the  angular 
limits  for  the  antenna  patterns.  OPCNST  is  used  to  initialize  various  constants  used  in  the  optical 
region  processing,  and  DCONST  performs  the  same  task  for  the  diffraction/troposcatter  region 
calculations. 

Once  the  initialization  phase  is  completed,  the  next  step  is  to  determine  if  the  input  range  is  in 
the  optical,  intermediate,  or  diffraction  region.  Subroutine  OPLIMIT  is  used  to  calculate  the  optical 
region  maximum  range  and  the  value  of  Ld  which  is  used  to  determine  Fin  the  diffraction  and 
intermediate  regions.  OPLIMIT  makes  use  of  the  ANTPAT,  OPFFAC,  R I  ITER,  REF,  and  RUFF 
subroutines.  OPFFAC  is  used  to  calculate  most  of  the  terms  used  in  antenna  pattern  factors, /(t,), 
REF  returns  R  and  4>,  and  RUFF  returns  the  surface  roughness  coefficient.  R I  ITER  is  used  to 
determine  the  reflection  point  range.  If  the  input  range  is  less  than  the  optical  region  limiting  range, 
subroutine  OPT1CF  is  called  to  determine  the  exact  value  of  the  pattern  propagation  factor,  F,  to  be 
returned.  This  requires  a  solution  to  the  cubic  equation  (Eq.  9)  to  determine  the  reflection  point 
range.  OPTICF  uses  a  Newton  iteration  technique  to  determine  the  reflection  point  range  and  then 
utilizes  OPFFAC,  ANTPAT,  REF,  and  RUFF  to  return  the  value  of  the  pattern  propagation 
factor,  F. 

If  the  input  range  is  greater  than  the  optical  region  maximum  range,  it  is  compared  to  rd,  the 
minimum  range  where  the  diffraction  region  calculations  are  valid.  If  the  range  value  is  greater  than 
rd.  then  the  DIFF  subroutine  is  used  to  determine  the  value  of  Ffor  this  region  D1FF  obtains  the 
appropriate  F  value  by  calling  the  HGAIN  and  TROPO  subroutines.  HG  AIN  returns  the  value  of 
the  height-gain  function  for  both  evaporation  and  surface-based  ducts  from  elevated  layers  and 
TROPO  returns  Ffor  the  troposcatter  region.  DIFF  calculates  the  value  of  Ffor  either  a  surface- 
based  duct  or  an  evaporation  duct.  If  a  surface-based  duct  containing  the  transmitter  exists,  it  is 
assumed  to  be  the  dominant  propagation  mechanism  and  Fis  calculated  by  using  Eq.  83  and  Eq.  3. 
If  no  surface-based  duct  exists,  then  F  is  determined  by  using  Eq.  57  and  Eq.  3.  The  diffraction 
region  Fis  then  compared  to  the  troposcatter  region  F.  If  the  troposcatter  F  value  is  within  18  dB. 
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the  fields  are  added  and  the  resulting  Fis  returned.  If  the  troposcatter  F  is  18  dB  less  than  the 
diffraction  F  value,  then  troposcatter  is  the  dominant  propagation  mechanism  and  the  troposcatter 
Fis  returned. 

If  the  input  range  is  less  than  rd  and  greater  than  the  optical  region  maximum  range,  then  the 
range  is  in  the  intermediate  region  and  a  linear  interpolation  of  the  propagation  factor  in  decibels 
versus  range  is  performed  to  obtain  F.  This  is  accomplished  by  calling  the  DIFF  subroutine  to 
obtain  the  value  of  Fat  rd.  Then,  a  linear  interpolation  between  the  value  of  Fat  the  optical  region 
maximum  range,  obtained  from  OPLIM1T,  and  Fat  rd  is  performed  to  obtain  the  value  of  Fat  the 
specified  range. 

FFACTR  is  structured  to  return  a  single  value  of  Ffor  a  given  set  of  inputs.  This  structure  is 
the  most  efficient  one  if  the  input  data  are  going  to  vary  for  each  case  in  some  arbitrary  fashion. 

This  is  not  the  case,  for  example,  in  the  demonstration  program,  MAIN,  where  only  the  range 
varies.  Since  the  environment  and  the  rest  of  the  EM  system  parameters  remain  constant  for  each 
call  to  FFACTR,  there  is  a  fair  amount  of  redundancy  in  the  calculation  of  the  antenna  pattern 
constants,  the  effective  earth  radius,  and  the  optical  and  diffraction  region  constants.  If  FFACTR 
is  always  going  to  be  used  in  such  a  fashion,  it  would  be  best  to  remove  these  calculations  from 
FFACTR  and  place  them  in  the  calling  routine.  For  example,  if  the  environmental  parameters 
are  constant  and  the  EM  system  height,  //,,  is  a  constant,  then  the  effective  earth  radius  is  also  a 
constant.  In  this  case,  the  DUCTS,  GETK,  INSRT,  MPROF,  and  PUSH  subroutines  should  be 
placed  in  the  calling  routine.  If  the  EM  system  is  fixed,  then  the  antenna  constants  subroutine, 
ANTPAR,  should  be  in  the  calling  routine.  If  the  environment  is  a  constant  tnd  the  only  indepen¬ 
dent  EM  system  variable  is  range,  as  in  the  MAIN  routine  supplied,  then  the  DCONST,  OPCNST, 
and  OPLIM1T  subroutines  should  also  be  placed  in  the  calling  routine.  This  would  mean  that 
ANTPAR.  DCONST.  DUCTS.  GETK.  INSRT,  MPROF.  OPCNST.  and  PUSH  subroutines 
would  only  be  called  once  for  each  environment  and  geometry  instead  of  the  multiple  calls  that 
occur  for  each  range  as  the  FFACTR  routine  is  presently  structured.  If  the  environment  is  to  remain 
constant,  but  both  the  radar  target/receiver  height  and  the  range  can  vary  in  an  arbitrary  fashion, 
then  the  DCONST,  OPCNST,  and  OPLIM1T  subroutines  would  have  to  remain  in  FFACTR,  but 
the  others  could  be  removed  to  the  calling  routine. 
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4.0  TEST  CASES 


A  number  of  EM  system  and  environmental  inputs  are  required  to  determine  the  propagation 
factor.  Table  4  !  lists  EM  system  parameters  for  five  test  case  systems,  Sysl  through  Sys5,  which  are 
used  to  verify  the  proper  operation  of  the  FFACTR  program.  The  input  parameters  that  are  listed 
correspond  to  the  variable  names  of  Section  3.0,  except  for  the  polarization  and  antenna  type 
entries.  Range  is  not  listed  as  an  EM  system  input  in  Table  4-1,  since  the  demonstration  program. 
MAIN,  supplies  a  constantly  decreasing  range  that  varies  from  100  km  to  (0  km  in  10-km  incre¬ 
ments  for  each  EM  system  and  environmental  data  set  selected.  Three  different  environmental  test 
case  conditions,  Envl  through  Env3,  are  listed  in  Table  4-2.  The  environment  of  Envl  corresponds 
to  a  standard  atmosphere  A/-unit  profile,  a  0-meter  evaporation  duct  height,  and  10  knots  of  wind. 
Env2  uses  the  same  Af-unit  profile  as  Envl,  but  the  wind  speed  is  0  knots,  and  a  10-meter  evapora¬ 
tion  duct  is  present.  Env3  has  an  Af-unit  profile  that  contains  a  100-meter  surface-based  duct, 
0-meter  evaporation  duct,  and  5  knots  of  wind.  The  environmental  data  sets  also  list  the  parameters 
using  the  variable  names  of  Section  3.0.  Each  of  the  EM  system  test  cases  uses  each  of  the 
environments. 


Table  4-1.  EM  system  test  set  input  data. 


Parameter 

Sysl 

Sys2 

Sys3 

Sys4 

Sys5 

/.  MHz 

100.0 

300.0 

5000.0 

9600.0 

15,000.0 

//,.  m 

20.0 

300 

200 

1.0 

50.0 

Hr  m 

30000.0 

2000.0 

200 

20.0 

120.0 

Polarization3 

H 

V 

c 

H 

H 

Antenna  Typeb 

S 

c 

H 

S 

O 

B  W,  deg 

10.0 

10.0 

2.0 

2.0 

n/a 

M0.  deg 

0.0 

1.0 

0.0 

1.0 

n/a 

aH  =  horizontal  bS-sin(x)'x 

V  -  vertical  C  -  esc2 

C  =  circular  H  =  generic  height-finder 


O  -  omnidirectional 
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Table  4-2.  Environmental  test  set  input  data. 


Parameter 

Envl 

Env2 

Env3 

6,  m 

0.0 

10.0 

0.0 

W,  kt 

10.0 

0.0 

5.0 

H,,  M,.( m,  M) 

(0.0,  339.0) 

(0.0,  339.0) 

(0.0.  350.0) 

H;,  M2,  (m,  M) 

(1000.0,  457.0) 

(1000.0, 457.0) 

(270.0,  381.9) 

H,,  M},  (tn,  M) 

(10000.0,1519.0) 

(10000.0,1519.0) 

(300.0,  340.0) 

H,,  M4,  (m,  M) 

n/a 

n/a 

(1000.0,  422.6) 

Hs,  My  <m,  M) 

n/a 

n/a 

(10000.0,1484.6) 

Tables  4-3  through  4-5  list  the  expected  output  data  for  the  different  environmental  test  cases. 
The  outputs,  in  decibels,  are  listed  to  the  nearest  0.1  dB,  and  the  FFACTR  program  is  considered  to 
be  operating  correctly  if  the  output  is  within  0. 1  dB  of  the  value  listed  in  the  appropriate  table. 


Table  4-3  Output  data  for  environment  I. 


Range  (km) 

fldB) 

Sysl 

Sys2 

Sys3 

Sys4 

Sys5 

100.0 

-37.2 

-0.6 

-56.1 

-62.5 

-58.0 

90.0 

-25.2 

-5.7 

-55.5 

-61.8 

-50.2 

80  0 

39.4 

-1.6 

-54.9 

-6i.2 

-28.1 

70.0 

-25. 8 

2.4 

-54.1 

-60.5 

-4.7 

60.0 

-25.3 

1.4 

-49.9 

-59.8 

0.8 

50.0 

-44.6 

-2.1 

-37.1 

-59.0 

2.7 

40.0 

-47.4 

1.0 

-23.1 

-53.7 

2.7 

30.0 

-26.2 

1.0 

-9.1 

-38.2 

0.3 

20.0 

-49.8 

1.0 

3.5 

-22.6 

-11.7 

10.0 

-43.2 

-0.1 

-1.6 

-7.7 

-4.0 

I 

i 
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Table  4-4.  Output  data  for  environment  2. 


HdB) 

Range  (km) 

Sy*l 

Syt2 

Sy*3 

Sy*4 

Sys5 

100.0 

-37.1 

-0.6 

-39.3 

-25.5 

-11.3 

90.0 

-25.1 

-5.7 

-33.7 

-22.9 

-11.8 

80.0 

-39.3 

-1.6 

-28.1 

-20.4 

-8.1 

70.0 

-25.7 

2.5 

-22.5 

-18.0 

-l.l 

60.0 

-25.1 

1.4 

-17.0 

-15.6 

0.8 

50.0 

-44.6 

-2.1 

-1 1.9 

-13.4 

2.7 

40.0 

-47.9 

1.0 

-7.3 

-11.3 

2.8 

30.0 

-25.9 

1.0 

-2.8 

-II. 1 

0.6 

20.0 

-61.7 

1.0 

3.5 

-12.2 

-18.9 

10.0 

-44.3 

-0.1 

-1.6 

-7.7 

-24.7 

Table  4-5.  Output  data  for  environment  3. 

F\  dB) 

Range  (km) 

Syil 

Syi2 

Sy»3 

Sy»4 

Sy»5 

100.0 

-41.7 

-5.1 

2.9 

0.1 

5.6 

90.0 

-24,8 

-4,0 

2,9 

0.1 

1.8 

80.0 

-44.1 

-0.9 

2.9 

0.1 

1.5 

70.0 

-26.0 

2.8 

2.9 

0.1 

-4,9 

60.0 

-25.4 

0.2 

2.1 

0.1 

0.5 

50.0 

-50.5 

-1.2 

1.3 

0.1 

5.0 

40.0 

-55.4 

0.5 

0.5 

-3.2 

4.6 

30.0 

-25.7 

0.8 

-0.3 

-7,9 

5.3 

20.0 

-54.0 

0.9 

5.1 

-12.7 

5.6 

10.0 

-43.6 

-0.1 

1.6 

-6.6 

0.6 
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Appendix 

STANDARD  PROPAGATION  MODEL  PROGRAM  LISTING 


c 

C  'k'k-k'k'k'klc'k'k-k'ftic'klc'k'k'k'k'k'k'k'k'k'kir'kitic  MAIN  PROGRAM  ********  it  Kirk  *  it  itit'kit'k'k'k'kirk'k'k 


C 


c  INPUTS: 
c 

c  EM  SYSTEM: 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


VARIABLE 

f  req 

ht 

hr 

polar 

antype 


bwidth 

elevat 


r 


VARIABLE  DESCRIPTION  (VALID  RANGE,  UNITS) 

SYSTEM  FREQUENCY  (100  -  20000  MHz) 

TRANSMITTER  ANTENNA  HEIGHT  (1  -  100  m) 
RECEIVER/TARGET  HEIGHT  (1  -  30000  m) 

ANTENNA  POLARIZATION  (HORIZONTAL  -  "H" , 
VERTICAL  -  "V",  CIRCULAR  -  "C") 
ANTENNA  TYPE  (OMNIDIRECTIONAL  -  "0" , 

SIN(X)/X  -  "S".  COSECANT- SQUARED  -  "C" . 

GENERIC  HT-FINDER  -  "H") 
ANTENNA  BEAM  WIDTH  (.5  -  45  DEG) 

ANTENNA  ELEVATION  ANC (-10  -  +10  DEG) 

(0  DEGREES  IS  HORIZONTAL,  NORMAL  POINTING 
ANGLE  FOR  SHIPBOARD  RADAR  SYSTEMS) 

RANGE  AT  WHICH  F- FACTOR  IS  DESIRED  (1- 1000km) 


c 


c  ENVIRONMENTAL: 


c 

delta 

c 

height(i) 

c 

munits(i) 

c 

wind 

c 

EVAPORATION  DUCT  HEIGHT  (0  -  40  m) 
HEIGHT  ARRAY  IN  METERS  -  UP  TO  30  ELEMENTS 
M-UNIT  ARRAY  CORRESPONDING  TO  HEIGHT  ARRAY 
WIND  SPEED  (0  -  50  KNOTS) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  OUTPUT: 

ff  20*L0G10(  PATTERN  PROPAGATION  FACTOR  )  dB 

( f f  values  chat  are  positive  indicate  a 
signal  level  above  the  free -space  field 
at  that  range.  Negative  values  indicate 
signal  levels  below  the  free-space  field. ] 

The  following  program  is  a  demonstration  driver  for  the  FFACTR 
(sub) routine .  This  program  is  included  to  show  possible  uses  for 
the  FFACTR  program.  The  FFACTR  program  returns  a  value  (in  dB) 
of  20*LOG(F;  where  F  is  the  pattern  propagation  factor.  F  is 
defined  as  the  ratio  of  the  actual  field  at  some  point,  to 
the  free-space  field  at  that  same  point.  (The  free-space  field 
is  determined  for  an  isotropic  antenna.)  Because  FFACTR  can  be 
called  in  any  arbitrary  fashion  It  is  not  necessarily  the  most 
efficient  structure  for  producing  a  product  such  as  a  loss-versus- 
range  plot.  If  only  the  range  is  to  be  varied,  with  constant 
terminal  heights,  then  the  ANTPAR,  DCONST,  DUCTS,  GETK ,  INSRT, 
MPR0F,  0PCNST  and  PUSH  subroutines  should  be  moved  to  the  calling 
routine  so  they  aren't  called  for  every  new  range  point. 


START  DEMO  PROGRAM 


include  ' envsys . common' 
c 
c 

real*4  ff,  r,  rloss 
integer*2  ZW 
c 

c  Enter  the  environmental  parameters, 

c 

call  envinp(delta ,  height,  Munits,  wind,  nmax) 
c 

c  Enter  the  EM  system  parameters. 
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c 

call  ays lnp ( f req , ht , hr , polar , antype , bvldch , elevat ) 
c 

fsterm  -  32.45  +  20.0  ★  ALOGlO(freq) 
dr  -  10.0 
r  -  100.0  +  dr 
c 

c  Call  FFACTR  for  the  EM  system  and  environmental  parameters 

c  entered  above.  Calculate  the  propagation  factor  for  several 

c  ranges  in  the  loop  below.  Print  the  output  values  of  prop- 

c  agation  factor  and  propagation  loss  for  each  range, 

c 

ZU  -  6 

DO  i  -  1,  10 
r  -  r  ■  dr 
call  ffactr(r,  FF) 

rloss  -  fsterm  +  20.0  *  AL0G10(r)  -  ff 
vrite(ZW,1000)  r.rloss.ff 

1000  Format ( 'Range  -  \f5.1,'  ka.  Loss  -  \f6.2,'  dB\'  F  - 
1  f6 . 1 , '  dB' ) 

END  DO 

END  DEMO  PROGRAM 

c 

END 
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Subroutine  ffactr 


FFACTR  returns  the  value  of  the  pattern  propagation  factor,  F,  in  dB 
for  specified  range  ,  EM  system  parameters  and  environmental  para* 
meters. 


Variables:  Description: 


alphac 

antbwr 
antelr 
antfac 
antype 

bwldth 
delta 
deltaf 

df  fac 
difac 
elevat 
elraaxr 
ff 

freq 
frsubd 
fzr 
fzt 

height 

hdif 

ht 
hr 
hi 
h2 

lvlant 
Munits 
nmax 
opmaxd 
opmaxf 
patrfac 
pd 

polar 

psi 
r 

rsdfac 
rsubd 
c 

c  sbdht 

c  theta 

c 

c  wind 

c 


Critical  angle  -  1st  angle  not  trapped  in  surface- 
based  duct. 

Antenna  vertical  beamvidth  in  radians . 

Antenna  elevation  angle  in  radians. 

Antenna  pattern  constant. 

Antenna  type:  0  -  omnidirectional,  S  -  sin(x)/x, 

C  -  cosecant-squared,  H  -  height- finder. 

Antenna  vertical  beam  width  in  degrees. 

Evaporation  duct  height  in  meters. 

Variable  used  in  linear  Interpolation  of  F  in  the 
intermediate  region. 

Diffraction  field  constant,  dB. 

Diffraction  field  constant,  dB. 

Antenna  elevation  angle  in  degrees. 

Maximum  elevation  angle  in  main  beam  of  antenna,  rad. 
Pattern  propagation  factor,  F,  in  dB. 

EM  system  operating  frequency  in  MHz. 

Pattern  propagation  factor  at  rsubd. 

Evaporation  duct  height-gain  at  hr,  dB. 

Evaporation  duct  height-gain  at  ht,  dB. 

Array  containing  environmental  input  height  values 
corresponding  to  the  Munit  array. 

Height  difference  between  receiver/target  and 
transmitter  heights  in  km. 

Transmitter  height  in  m. 

Receiver/target  height  in  m. 

Lower  height  of  hr,  ht,  in  m 
Higher  height  of  hr,  hr,  in  ra. 

Transmitter  height  level  in  hmrs  and  dMdh  arrays. 

Array  containing  environmental  input  H-unit  values. 
Integer  number  of  layers  in  Munits  and  height  arrays. 
Maximum  range  in  the  optical  interference  region,  km. 

F  at  opmaxd. 

Antenna  pattern  constant. 

Path-difference  between  direct  and  sea-reflected  rays. 
Antenna  polarization:  H  -  horizontal,  V  -  vertical, 

C  »  circular. 

Crazing  angle  in  radians. 

Range  in  km. 

Constant  used  to  calculate  rsubd. 

Minimum  range  where  diffraction  field  solutions  are 
applicable,  km. 

Surface-based  duct  height,  m. 

Total  phase  difference  between  direct  and  sea- 

reflected  rays  including  phase  lag  due  to  reflection. 
Wind  speed  in  kts . 


SUBROUTINE  ffactr(r,  FF) 
c 
c 

real*4  deltaf,  ff,  fzr,  fzt,  opmaxd,  opmaxf,  frsubd,  r 
real*4  dMdh(32),  hmrs(32) 
integer*2  lvlant,  ntot 
c 
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include  ' ffac . common' 
include  ' envsys . common' 
c 

c  Call  mprof  to  insert  a  profile  level  at  Ht  and  determine  if 

c  any  surface-based  ducts  are  present.  If  a  surface-based  duct 

c  is  present  calculate  critical  angle,  alphac. 

c 

call  mprof (height,  Munits,  ht,  NMAX,  ALPHAC,  DMDN,  HMRS, 

1  SBDHT,  NTOT) 

c 

c  Call  getk  to  determine  the  effective  earth  radius  factor,  rk. 

c 

call  getk(alphac,  dMdh,  hmrs,  ntot,  ht,  RK) 
c 

c  Define  hi,  h2  for  opticf  subroutine.  These  are  swapped  for 

c  ht>hr  because  the  iteration  loop  for  rl  in  opticf  works  most 

c  efficiently  when  the  lowest  height  is  the  transmitter  height, 

c 

IF  (ht  .GT.  hr)  THEN 
hi  -  hr 
h2  -  ht 
ELSE 
hi  -  ht 
h2  -  hr 
END  IF 

c  Define  optical  region  constants, 

call  opcnst 

hdif  -  (hr  -  ht)  *  1.0e-3 
c  Initialize  antenna  parameters. 

call  antpar ( antype , bwidth , e levat , ANTBWR , ANTELR , ANTFAC , 

1  ELMAXR , PATRFAC ) 

c  Define  diffraction/troposcatter  region  constants, 

call  dconst 
call  hgain(hr ,  FZR) 

IF  (sbdht  .EQ.  0.0)  THEN 
call  hgain(ht ,  FZT) 
dffac  -  dffac  -  fzt 
END  IF 

difac  -  dffac  -  fzr 

rsubd  -  3.572  *  (SQRT(rkmin  *  ht)  +  SQRT(rkmin  *  hr))  +  rsdfac 
c  Determine  maximum  range  and  f-factor  in  optical  region, 

call  oplimit(OPMAXD,OPMAXF) 

IF  (r  .GE.  rsubd)  THEN 

c  Calculate  loss  for  range  in  diffraction/troposcatter  region, 

call  diff (r ,  FF) 

ELSE 

IF  (r  .GT.  opmaxd)  THEN 

Ranee  is  in  intermediate  region  •  use  linear  interpolation 
on  log  of  the  f-factor. 
call  diff( rsubd, FRSUBD) 

deltaf  -  (r  -  opmaxd)  *  (opmaxf  -  frsubd)  /  (opmaxd- rsubd) 
ff  -  opmaxf  +  deltaf 
ELSE 

c  Range  is  in  the  optical  interference  region. 

IF  (r  ,LE.  opmaxd)  call  opticf (polar , r , PD, PSI , THETA, FF) 

END  IF 
END  IF 
ff  -  -  ff 
RETURN 
END 
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C 

c  'envsys. common'  include  file 

c 

c  EM  system  parameter  common  blocks 

c 

common  /  emsystem  /  freq,  hr,  ht 
common  /  emsystem  /  polar,  antype,  bwidth,  elevat 
c 

c  Environmental  parameter  common  blocks 

c 

common  /  enviro  /  delta,  height,  Munits,  nmax,  wind 
c 
c 

real*4  delta,  height(30),  Munits(30),  wind 
real*4  freq,  ht,  nr,  bwidth,  elevat 
character*l  antype,  polar 
integer*2  nmax 
c 


c 

c 

c 


c 


' ffac . common'  include  file 

common  /  comffactr  /  ae ,  ae2 ,  aeth,  alpha,  alphac,  antbwr 
common  /  comffactr  /  antelr,  antfac,  atten 
common  /  comffactr  /cl,  c2  ,c3  ,  c4  ,c5,  c6,  cl 
common  /  comffactr  /  del,  dffac,  difac,  elmaxr,  exloss 
common  /  comffactr  /  fsterm,  hbar,  hbfreq,  hdif,  hmin,  hi 
common  /  comffactr  /  h2 ,  horznl,  patd,  patrfac,  rk,  rkmin 
common  /  comffactr  /  rnimag,  rnreal,  rsdfac,  rsubd,  sbdht 
common  /  comffactr  /  thefac,  twoae ,  zfac,  zmax 

real*4  ae,  ae2  aeth,  alpha,  alphac,  antbwr,  antelr,  antfac, 
atten,  cl,  c2,  c3,  c4,  c5,  c6,  cl,  del,  dffac,  difac, 
elmaxr,  exloss,  fsterm,  hbar,  hbfreq,  hdif,  hmin, 
horznl,  hi,  h2 ,  patd:  patrfac,  rk,  rkmin,  rnimag, 
rnreal,  rsdfac,  rsubd,  sbdht,  thefac,  twoae,  zfac,  zmax 
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c 

c 

c 

c 

c 

c 


Subroutine  antpar 

ANTPAR  is  used  to  intitialize  antenna 
calculating  antenna  pattern  factors. 


Description: 


parameters  for  use  in 


c 

Variable : 

c 

c 

Antbwr 

c 

Antelr 

c 

Antfac 

c 

Antype 

c 

c 

c 

c 

Bwidth 

c 

Elevat 

c 

Elmaxr 

c 

Patrfac 

Antenna  bean  width  in  radians . 

Antenna  elevation  angle  in  radians. 

Antenna  pattern  constant. 

Antenna  type:  0  -  omnidirectional 
S  -  Sin(x)/x 
C  •  Cosecant-squared 
H  -  generic  Height- finder 
Antenna  beam  width,  degrees. 

Antenna  elevation  angle,  degrees. 

Maximum  angle  in  main  beam  of  antenna,  radians. 
Pattern  factor  constant  for  Sin(x)/x  antennas, 
used  to  calculate  Elmaxr  for  Sin(x)/x  antennas. 

SUBROUTINE  antpar ( antype , bwidth , e leva t , ANTBWR , ANTELR , ANTFAC , 
1  ELMAXR, PATRFAC) 


real*4  antbwr,  antelr,  antfac,  amax,  bwidth,  elmaxr, 

1  elevat,  pi,  patfac,  patrfac 

character*l  antype 

PI  -  3.14159 

c  Convert  beam  width  and  elevation  angle  to  radians, 

antbwr  -  1 . 745e- 2*bwidth 
antelr  -  1 . 745e - 2*elevat 
elmaxr  -  1.047 
IF  (antype  .NE.  "0")  THEN 
IF  (antype  .EQ.  "C")  THEN 

c  Cosecant-squared  antenna  pattern  constants, 

elmaxr  -  antelr  +  .78525 
antfac  -  SIN(antbwr) 

ELSE 

IF  ((antype  .EQ.  "S") .OR. (antype  .EQ.  ”H" ) )THEN 
c  Sin(x)/x  and  height-finder  antenna  pattern  constants, 

antfac  -  1 . 39157/SIN(antbwr/2 .0) 
amax  -  Pl/antfac 

patrfac  -  -ATAN(amax/SQRT(l . 0  -  amax*amax)) 

IF  (antype  .EQ.  "S")  elmaxr  -  antelr  -  patrfac 
END  IF 
END  IF 
END  IF 
RETURN 
END 
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Subroutine  antpat 


ANTPAT  returns  the  antenna  pattern  factor  for  a  given  angle 
and  antenna  type . 


Variable : 
alpha 
antbwr 
antelr 
antfac 
angle 
antype 


patfac 

patrfac 


Description: 

Direct  r'y  launch  angle,  radians. 

Antenna  beam  width  in  radians. 

Antenna  elevation  angle  in  radians. 

Pattern  constant. 

The  angle  for  which  the  pattern  factor  is  desired. 
Antenna  pattern  type:  0  -  omnidirectional 

S  -  sin(x)/x 
C  -  cosecant* squared 
H  -  generic  height- finder 
The  antenna  pattern  factor  for  the  given  angle. 
Pattern  constant. 


SUBROUTINE  antpat(antype .alpha, antbwr .antelr, antfac .patrfac , 
1  angle, PATFAC) 


real*4  alpha,  alphaO,  antbwr,  antelr,  antfac,  angle,  apat, 

1  patfac,  patrfac,  ufac 

character*l  antype 

c 

patfac  -  1.0 

IF  (antype  .NE.  "0")  THEN 

c  Antenna  types  other  than  omni  require  calculation. 

IF  ((antype  .EQ.  "H") .AND. (alpha  .GT.  antelr))  THEN 
alphaO  -  alpha 
ELSE 

alphaO  -  antelr 
END  IF 

apat  -  angle  -  alphaO 
IF  (antype  .EQ.  "C")  THEN 
c  Cosecant -squared  antenna  type. 

patfac  -  AMIN1U.0,  AMAX1(0 . 03 ,  1.0  +  apat/antbwr) ) 

IF  (apat. GT. antbwr)  patfac  -  SIN(antbwr)/SIN(ABS(apat) ) 
ELSE 

c  SIN(X)/X  antenna  type. 

IF  (apat  NE.  0.0)  THEN 

IF  ((angle  .LE.  alphaO  +  patrfac)  .OR. 

1  (angle  .GE.  alphaO  -  patrfac))  THEN 

c  Antenna  pattern  is  limited  to  main  lobe  only, 

patfac  -  0.03 
ELSE 

c  Sin(x)/x  calculation. 

ufac  -  antfac*SIN(apat) 

patfac  -  AMINK1.0,  AMAX1(0.03,  SIN(ufac)/ufac) ) 
END  IF 
END  IF 
END  IF 
END  IF 
RETURN 
END 
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c 

c  Subroutine  dconst 
c 

c  DCONST  initializes  variables  for  the  diffraction  and  troposcatter 
c  region  routines, 
c 

c  Variable: 
c  arg 

c  atten 

c  cl  -  c7 

c  del 

c  delta 

c  dffac 

c  fmax 

c  freq 

c  fsterm 

c  gamma 

c  hmin 

c  rfac 

c  rk 

c  rkmin 

c 

c  rsdfac 

c  sbdht 

c  zfac 

c  zraax 

c 
c 
c 

SUBROUTINE  dconst 


real*4  arg,  fmax,  gamma,  rfac,  slope 
c 

include  ' ffac . common' 
include  ' envsys . common' 
c 
c 

c  Diffraction  region  constants. 

c 

IF  (sbdht  .CT.  0.0)  THEN 
c 

c  Surface -based  duct  model, 

c 

del  -  0.0 
hmin  -  1.0 
atten  -  0.0 
dffac  -  fsterm 
c 

ELSE 

c 

c  The  following  terms  are  for  NOSC  evap  duct  model, 

c 

rfac  -  0.04705  *  freq**( 1 . /3 . ) 
zfac  -  0.002214  *  freq**( 2 . /3 . ) 
hmin  -  1.0 

del  -  AMINl(delta  *  zfac,  23.3) 
c 

c  Constants  for  scaled  evap.  duct  heights  >—  10.25  meters. 

cl  -  -0.1189  *  del  +  5.5495 
c3  -  3  .  / 2  . 

c2  -  1.3291  *  SIN(0 . 218  *  (del- 10 . 0)**0 . 77)  +  0 . 2171*ALOC(del ) 
c2  -  c2  *  4 . /2**( -c3) 


Description: 

Evaporation  duct  model  temporary  variable. 

Diffraction  region  attenuation  rate  in  dB/km. 
Evaporation  duct  constants  for  height-gain  function. 
Scaled  evaporation  duct  height  (delta  *  zfac). 
Evaporation  duct  height,  m. 

Diffraction  field  constant  in  dB. 

Evaporation  duct  model  temporary  variable. 

EM  system  frequency  in  MHz. 

Free -space  loss  term,  dB. 

Evaporation  duct  excitation  factor  in  dB. 

Minimum  allowable  height,  m. 

Evaporation  duct  range  scale  factor. 

Effective  earth  radius  factor. 

Minimum  rk  used  for  calculation  of  the  diffraction 
region  minimum  range,  rsubd. 

Constant  used  for  rsubd  calculation. 

Surface -based  duct  height,  m. 

Evaporation  duct  height  scale  factor. 

Evaporation  duct  height  variable.  Height  where  the 
two  different  equations  for  the  height-gain  factors 
must  be  equal  (del  >—  10.25  meters). 
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c4  -  87.0  -  SQRT(313 . 29  -  (del  -  25.3)**2) 
zmax  -  4.0  *  EXP(-0. 31*(del  -  10.0))  +  6.0 
arg  -  c2  *  zmax**c3 

slope  -  4 . 72  *  cl  *  c2  *  c3  *  SQRT(zmax)  /  TAN(arg) 
c7  -  49.4  *  EXP( -0 . 1699* (del  -  10.0))  +  30.0 
fmax  -  cl  *  AL0G(SIN(arg) )  +  c4  -  c7 
c6  -  (zmax/4.72)  *  slope  /  fmax 
c5  -  fmax  /  zmax**c6 
ELSE 

Constants  for  scaled  evap.  duct  heights  <-  10.25  meters. 

c2  -  SQRT(40623 .61  -  (del  +  4.4961)**2)  -  201.0128 
cl  -  (-2.2  *  EXP(-0.244*del)  +  17 .0)*4.72**(-c2) 
c4  -  SQRT( 14301 . 2  -  (del  +  5.32545)**2)  -  119.569 
c3  -  (-33.9  *  EXP( -0 . 5170001*del)  -  3 . 0)*4 . 72**( -c4) 
c5  -  41.0  *  EXP( -0 . 41*del)  +  61.0 
END  IF 

atten  -  92.516  -  SQRT(8608 . 7593  -  (del  -  20.2663)**2) 

IF  (atten  .LT.  0.0009)  atten  -  0.0009 
atten  -  atten  *  rfac 

IF  (del  .LE.  3.8)  gamma  -  216.7  +  del  *  1.5526 
IF  (del  .CT.  3.8)  gamma  -  222.6  -  (del  -  3.8)  *  1.1771 
dffac  -  51.1  +  gamma  +  10.0  *  ALOGlO(rfac) 

END  IF 

Constants  used  to  calculate  rsubd,  the  minimum 

range  at  which  the  diffraction  field  solutions  are  valid. 

rkmin  -  AMAXl(rk,  1.3333) 

rsdfac  -  230.2  *  (rkmin**2  /  freq)**(l. 0/3.0) 

RETURN 

END 


Subroutine  diff 

Subroutine  DIFF  returns  the  diffraction  field  propagation  factor 
as  a  function  of  range. 


VARIABLES : 
atten 
delta 
df loss 
dloss 
dif 
difac 
diffe 
exloss 
r 

tloss 


DESCRIPTION: 

NOSC  model  attenuation  rate  in  dB/km 
Evaporation  duct  height  in  meters 
-  20*LOG(F) ,  where  F  is  the  propagation  factor 
Diffraction  field  strength  in  dB 
Temporary  variable 

NOSC  evaporation  duct  model  constant 

NOSC  evaportion  duct  model  loss  in  dB 

Antenna  loss  for  lowest  angle  in  optical  region  (dB) 

range  in  km 

Troposcatter  loss  from  Tropo  Subroutine  in  dB 


c 

SUBROUTINE  diff(r.  DFLOSS) 
c 
c 
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real*4  dlf,  dfloss,  dloss,  dif fe ,  r,  tloss,  tlr 
c 

Include  ' ffac . common' 
include  ' envsys . common' 
c 

tlr  -  10 . 0*ALOG10(r) 

IF  (sbdht  .EQ.  0.0)  THEN 

c  Calculate  the  evaporation  duct  loss, 

dloss  -  difac  +  tlr  +  atten*r 
ELSE 

c  Calculate  the  surface  based  duct  loss, 

dloss  -  difac  +  2.0*tlr 
END  IF 

dloss  -  dloss  +  exloss 
c 

c  Calculate  troposcatter  loss  and  compare  to  dloss.  If  the 

c  difference  is  +/-  18  dB  add  the  two  fields  together. 

c 

call  tropo(r, tloss) 
dlf  -  dloss  -  tloss 
IF  (dif  .GE.  18.0)  THEN 
c  Troposcatter  field  dominates, 

dloss  -  tloss 

ELS El F  (dif  .GE.  -18.0)  THEN 

c  Add  troposcatter  and  diffractions  fields  togather. 

dloss  -  dloss  -  10.0*ALOG10(1.0  +  10.0**(dif/10.0) ) 

END  IF 

c 

c  -20*L0G(F)  -  actual  loss  -  free  space  loss 

c 

dfloss  -  dloss  -  fsterm  -  2.0*tlr 

RETURN 

END 


c 

c  Subroutine  ducts 

c 

c  DUCTS  builds  an  array  containing  the  top,  bottom,  and 
c  minimum  refractivity  of  all  the  major  ducts  in  the 
c  atmosphere  refractivity  profile. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Variable : 
dct 


lvls 

ndcts 


rht 


rmu 


Description: 

3,*  duct  parameters  array. 

1, n  bottom  of  duct  'n' ,  meters. 

2, n  top  of  duct  'n',  meters. 

3, n  minimum  refractivity  of  duct  'n',  M-units. 
Number  of  refractivity  level  in  rmu,  rhts. 

in:  the  maximum  number  of  ducts  allowed, 
out:  the  number  of  ducts  found. 

Duct  counter. 

Height  array,  meters. 

Modified  refractivity,  M-unlt  array,  elements 

correspond  to  like-number  elements  of  rht  array. 
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SUBROUTINE  ducts ( rmu , rht , lvis , DCT.NDCTS ) 
c 

real*4  dct,delu,delh,deltu,hbot,htop,rht(32) ,rmu(32) 
integer*2  lvls ,  ibot,  Lduct ,  lend,  iq,  itop.ndcts .nq 
dimension  dct(3,8) 

Locate  all  major  ducts 
nq-0 

iq-3*ndcts 

itop-lvls 

iend-ndcts 

ndcts-0 

DO  iduct-l,iend 

Look  for  top  of  next  duct 
1010  continue 

htop-rht( itop) 
if(itop.eq.l)  go  to  1060 
ibot-itop-1 

if (rmu(itop) . le . rmu(ibot) )  go  to  1020 
itop- itop- 1 
go  to  1010 

Look  for  bottom  of  the  duct 
1020  continue 

hbot-rht( ibot) 

if (rmu(ibot) . It . rmu( itop) )  go  to  1030 
if ( ibot . eq . 1)  go  to  1040 
ibot-ibot- 1 
go  to  1020 

Calculate  bottom  of  duct  using  linear  interpolation 
1030  continue 

dela-rmu(ibot+l) -rmu(ibot) 
delh«rht(ibot+l) -rht (ibot) 
deltu-rmu( itop) -rmu (ibot) 
if(delu.lt.O.Ol)  go  to  1040 
hbot-rht( ibot)  +  aeltu*delh/delu 
c 

c  Store  duct  parameters  in  array  dct 

1040  continue 

amu-rmu(itop) 

call  push(dct , iq , nq , amu) 

call  push(dct , iq ,nq ,htop) 
call  push(dct , iq ,nq ,hbot) 
ndcts-iduct 
itop-ibot 
END  DO 
c 

1060  continue 
RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 


c 


c 

c 


c 


c 


c 


c 


Subroutine  envfil 


ENVFIL  lists  the  available  environmental  files  and  allows  the 
user  to  select  one.  The  selected  environmental  file  is  read 
and  closed.  The  data  from  the  file  is  returned  to  the  calling 
routine . 


Variable : 
delta 
height 

levels 

Hunits 

wind 


Description: 

Evaporation  duct  height  in  m. 

Array  of  up  to  30  elements  containing  the  heights 
of  the  M-unit  profile. 

The  number  of  levels  in  the  height,  Munits  arrays. 
Array  of  up  to  30  elements  containing  the  M-unit 
values  of  the  upper-air  profile. 

Wind  speed  in  knots. 


SUBROUTINE  envfil (delta ,  height,  Munits,  wind,  levels) 


real*4  delta,  height(30),  Munits(30),  wind 
integer*2  levels,  ZR,  ZW 
character*12  filename 

Initialize  read,  write  channels 
ZR  -  5 
ZW  -  6 

write  (Z U,'("  Available  Environmental  Files:  ")') 

List  all  files  beginning  with  "E". 
call  system  ('Is  (EJ*  1>&2 '//char(0) ) 
write  (ZW, '(//, "Enter  input  file  name:  ",$)') 
read  (ZR,'(al2)')  filename 
open  (10,FILE-filename) 

Read  wind  speed  in  knots  and  evaporation  duct  height  in  ra. 
read  (10,  '(f4.1)')  delta 
read  (10,  '(f4.l)')  wind 

Read  the  number  of  levels  in  M-unit  profile, 
read  (10,  ' ( i2) ' )  levels 

Read  the  height  and  M-unit  profile  array  values. 

DO  i-1 ,  levels 

read  (10,  '(2fl0.1)')  height(i),  Munits(i) 

END  DO 

Close  environmental  file. 
close(10) 

RETURN 

END 


A-l  2 


o  o 


c 

c  Subroutine  envinp 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 


c 


Subroutine  ENVINP  prompts  the  user  to  enter  environmental  parameters 
and  returns.  Environments  can  be  entered  over  the  keyboard  or  from 
a  file.  If  the  environment  is  entered  over  the  keyboard  it  can  be 
saved  in  a  file  for  future  use. 


Variable : 
delta 
height 

levels 

Munits 

wind 


Description: 

Evaporation  duct  height  in  m. 

Array  of  up  to  30  elements  containing  the  heights 
of  the  M-unit  profile. 

The  number  of  levels  in  the  height,  Munits  arrays. 
Array  of  up  to  30  elements  containing  the  M-unit 
values  of  the  upper-air  profile. 

Wind  speed  in  knots. 


SUBROUTINE  envinp (DELTA,  HEIGHT,  MUNITS,  WIND,  LEVELS) 

real*4  delta,  height(30).  Munits(30),  wind 
character*20  A,  dummy,  filename 
integer*2  k,  kt,  levels,  ZW,  ZR 

Specify  the  read  (5)  and  write  (6)  channel  numbers. 

ZW  -  6 
ZR  -  5 

Initialize  environmental  parameters, 
wind  -  0.0 
delta  -  0.0 
levels  -  2 
DO  i  -  1,30 

height(i)  -  0.0 
Munlts(i)  -  0.0 
END  DO 

Enter  the  environmental  data  parameters. 
write(ZW, ' ("Enter  environmental  data  parameters.  You  may  enter")') 
write(ZW, ' ("up  to  30  layers  or  enter  data  from  a  file.  ")') 


c  Select  environmental  file. 

write(ZW, ' ("Enter  data  from  a  file?  (yes  or  no)  ",$)') 
read(zr, ' (A) ' )dummy 

IF  ((dummy(l : 1)  .eq.  ' y '). or . (dummy (1 : 1)  .eq.  'Y'))  THEN 
call  envfil(delta,  height,  Munits,  wind,  levels) 

ELSE 

write(ZW, ' ("Adjcent  layers  must  have  different  M-values  and")') 
write(ZW, ' ("at  least  two  layers  are  required.")') 

c 

height(l)  -  0.0 
Munlts(l)  -  0.0 
write(ZW, 1000) 

1000  format(/ ,' Enter  M-unit  Profile  -  (Height  in  meters,  M-units) ' 

1  /, 'Starting  height  is  at  surface  (0  meters)  ') 

c 

c  DO  loop  to  enter  profile  data  (Height  and  Munit  arrays). 


DO  i  -  1,  30 

100  write(zw,'("  Enter  height  in  meters  (or  end)  ",$)') 
read(zr , ' (A) ' )dummy 

IF  ( ( dummy (1:1)  .EQ.  ' e '). OR . (dummy ( 1 : 1)  .EQ.  'E'))  goto  200 
k  -  1 
kt  -  1 

DO  WHILE  ( (k.t  .eq.  11  .and.  (k  .le.  20)) 

IF  (dummy (k : k) . EQ . '  ')  dummy(k:k)  - 
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IF  (dummy (k : k) .EQ. kt  -  0 
k  -  k  +  1 
END  DO 

IF  (1  .gt.  1)  THEN 

read( dummy , ' (flO. 2) ' )height(l) 

IF  (height (1)  .LE.  height(i-l))  THEN 
write(zw,1010) 

1010  format ( 'Heights  must  increase,  re-enter  heifht  ') 

goto  100 
END  IF 
END  IF 
levels  -  i 

write(zw, ' ("  Enter  M*unit  value  at  level  ",$)') 
read(zr , ' (A) ' ) dummy 
k  -  1 
kt  -  1 

DO  WHILE((kt  .EQ.  1)  .AND.  (k  .LE.  20)) 

IF  (dummy(k:k) .EQ. '  ')  dummy(k:k)  - 
IF  (dummy(k:k) .EQ. ' . ' )  kt  -  0 
k  -  k  +  1 
END  DO 

read(dummy, ' (flO. 2) ' )Munits(i) 

IF  ( (i  ,NE.  1)  .AND.  (Munits(i)  .EQ.  Munits(i-l)))  THEN 
Hunits(i)  -  Munits(i)  +  0.1 
END  IF 
END  DO 

200  continue 

write(ZW, 1020) 
c 

1020  format( ' Enter  evaporation  duct  height  in  meters  (0  to  40)  ',$) 
read(ZR,*)  delta 
IF  (delta  .LT.  0.0)  delta  -  0.0 
IF  (delta  .GT.  40.0)  delta  -  40.0 
c 

write(ZW, 1030) 

1030  format( 'Enter  wind  speed  in  knots  (0  to  50)  '  ,$) 
read(ZR,*)  wind 
IF  (wind  .LT.  0.0)  wind  -  0.0 
IF  (wind  .GT.  50.0)  wind  -  50.0 
c 

write(ZW, ' ("Do  you  wish  to  store  this  environment  in  a  file?", 
1  "  (yes  or  no)  ",$)') 

read(zr , ' (A) ' )dummy 

IF  ( (duramy(l : 1)  .eq.  'y' ) .or . (dummy(l:l)  .eq  'Y'))  THEN 
write  (ZW,'(M  Current  Environmental  Files:  ")') 
call  system  ('Is  (EJ*  1>&2 '//char(0) ) 
write  (ZW,  1040) 

1040  format(" Enter  file  name  (First  letter  MUST  be  E)  ",$) 
read  (ZR.'(al2)')  filename 
open  (10 , FILE-filename) 
c 

c  Write  wind  speed  in  knots  and  evaporation  duct  height  in  m. 

write(10,  ' (f4 . 1) ' )  delta 
write(10,  ' (f4. 1) ' )  wind 

c  Write  the  numbers  of  levels  in  M-unit  profile. 

write(10.  '(12)')  levels 
DO  i— 1,  levels 

write(10,  ' (2f 10, 1) ' )  height(i),  Munits(i) 

END  DO 

c  close  file 

close( 10) 

END  IF 
c 

END  IF 
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c 


RETURN 

END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


c 


Subroutine  getk 

Subroutine  GETK  Is  used  to  determine  the  effective  earth  radius 
factor  k.  Getk  accomplishes  this  by  tracing  a  ray  from  the  trans¬ 
mitter  height  to  200  NM1  (370  km).  The  ray  launch  angle  is  0  deg. 
If  no  surface-based  duct  exists,  or  alphac,  the  critical  angle  if 
one  does. 


Variable : 
alphac 

aO 

al 

deld 

delh 

delM 

delmdh 

dMdh 

hlast 

hmrs 

ntot 

rdeld 

rmax 

rng 

rk 

xmtr 


Description: 

Critical  angle  necessary  to  escape  duct.  If  alphac 
-  0  then  no  surface -based  duct  exists. 

Initial  ray  launch  angle,  radians. 

Ray  angle  at  top  of  layer,  radians. 

Range  difference,  km. 

Height  difference,  meters. 

M-unit  difference. 

M-unit  gradient. 

M-unit  gradient  array. 

Height  at  370  km. 

Array  of  height  elements,  In  meters. 

Maximum  number  of  elements  in  hmrs  and  dMdh  arrays. 
Range  incremented  in  ray  trace. 

Maximum  range  for  ray  trace  -  370  km. 

Range ,  km . 

Effective  earth  radius  factor. 

Transmitter  height  in  meters. 


SUBROUTINE  getk( alphac,  dMdh,  hmrs,  ntot,  xmtr,  RK) 


real*4  alphac,  aO,  al,  deld,  delh,  delM,  delmdh,  dMdh(32) 
real*4  hlast,  hmrs(32),  rdeld,  rmax,  rng,  rk,  xmtr 
lnteger*2  ntot,  1 


rmax  -  370.0 
h  -  xmtr 
rng  -0.0 
aO  -  alphac 

Loop  to  trace  ray  through  the  atmospheric  layers. 
DO  1-2, ntot- 1 

delm  -  (hmrs(t+l)  -  h)*dMdh(i)*l .0E-3 
al  -  SQRT(a0*a0  +  2.0*delm) 
deld  -  (al  -a0)/dMdh(l) 
rdeld  -  rng  +  deld 
I F( rdeld  .CT.  rmax)  GOTO  1000 
aO  -  al 
h  -  hmrs(l+l) 
rng  -  rdeld 
END  DO 
1  -  ntot 
1000  continue 

Ray  trace  In  final  layer  to  range  rmax. 
ueld  -  rmax  -  rng 
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al  -  aO  +  dMdh(i)  *  deld 
delM  -  (al*al  -  a0*a0)*0.5 
delh  -  1000 . 0*delM/dMdh ( i) 
hlasc  -  hmrs(i)  +  delh 

c  Determine  the  equivalent  single -gradient  atmosphere  that 

c  would  be  required  to  trace  a  ray  launched  at  alphac  that 

c  would  arrive  at  height  -  hlast  at  a  range  of  370  km. 

delmdh  -  ( -alphac) *2 .0/rmax  +  2.0E-3*(hlast  -  xmtr)/(rmax*rmax) 
rk  -  1.0/(6371.0  *  deltndh) 

IF(rk  .GT.  5.0)  rk  -  5.0 
IF(rk  .LE.  0.5)  rk  -  0.50 

RETURN 

END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


Subroutine  gtheta 

GTHETA  calculates  optical  phase-lag  difference  angle  'theta' 
between  direct  and  sea- reflected  rays  using  the  reflection 
point  range  'rl' 


Variable : 
ae2 
hi 
h2 
hip 
h2p 
plr 


psl 

phi 

r 

rl 

r2 

rmag 

theta 


Description: 

Effective  earth  radius  *  2000. 

Height  of  transmitting  antenna,  m. 

Height  of  receiver/target,  m. 

Effective  height  of  hi,  m. 

Effective  height  of  h2,  m. 

Antenna  polarization:  H  -  horizontal 

V  -  vertical 
C  -  circular 

Crazing  angle  in  radians. 

Phase  lag  due  to  reflection,  radians. 

Total  ground  range,  km. 

Reflection  point  range,  (from  hi),  km. 
Reflection  point  range,  (from  h2),  km. 
Magnitude  of  the  reflection  coefficient. 
Total  phase  lag  between  direct  and  reflected 
rays  Including  phi. 


SUBROUTINE  gtheta (plr , r 1 , R , THETA , R2 , PSI , RMAG) 


real*4  hip,  h2p ,  psi,  phi,  r,  rl,  r2 ,  rmag,  theta 
character*!  plr 

include  ' ffac . common' 
inc lude  ' envsy s . common ' 

hip  -  hi  -  rl*rl/ac2 
psi  -  1.0e-3  *  hlp/rl 

IF  (psl  .GT.  0.03)  psi  -  ATAN ( 1 . Oe - 3  *  hlp/rl) 

Ray  trace  equation  used  to  determine  r2  based  on  psi. 
r2  -  (  SQRT(psi*psi  +  2.0e-3  *  h2/ae)  -  psl  )  *  ae 
r  -  rl  +  r2 
h2p  -  h2  -  r2*r2/ae2 
call  ref(plr , psi , RMAG, PHI) 
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c  Calculate  theta  -  Path* length  difference  +  phase  lag  due 

c  to  reflection  (phi) . 

theta  -  phi  +  thefac*hlp*h2p  /  r 

RETURN 

END 


c 

c 

c 

c 

c 

c 


c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


Subroutine  hgain 


HGAIN  returns  a  height-gain  factor  in  dB  for  a  specified  height. 


Variable: 
cl  -  c7 
del 
delta 
freq 
fzdb 
h 

hmin 

sbdht 

rfac 

zfac 

zmax 

zl 

z2 


Description: 

Constants  used  to  calculate  fzdb  for  evap.  ducts. 
Scaled  evaporation  duct  height. 

Evaporation  duct  height,  a. 

EM  system  frequency  in  MHz. 

Heignt-gain  factor  in  dB. 

The  height  for  whi  "  the  helght*galn  factor  is 
required,  m. 

Minimum  height. 

Surface -based  duct  height,  m. 

Evaporation  duct  range  scale  factor. 

Evaporation  duct  height  scale  factor. 

Height  above  which  different  eqn.s  are  used  for  height 

fain  calculations  for  del>lu.25m. 
ed  height  for  surface-based  ducts. 

Scaled  height  for  evaporation  duct  heights. 


SUBROUTINE  hgain  (h.  FZDB) 


real*4  fzdb,  h,  zl,  z2 


include  ' ffac. common' 

Include  'envsys. common' 

fzdb  -  0.0 

IF  (sbdht  .GT.  0.0)  THEN 

Calculate  surface-based  duct  height-gain  factor, 
zl  -  h  /  sbdht 

IF  ((Freq  .LE.  150 .0) . AND. (zl  .LT.  0.8))  THEN 
fzdb  -  -50.0  *  (zl  -  0. 5)**2 
END  IF 

IF  ((Freq  .LE.  150.0) .AND. (zl  .GE.  0.8))  THEN 
fzdb  -  1.14  *  zl**( -6.26)  -  10.0 
END  IF 

IF  ((Freq  .GT.  150.0) .AND. (zl  .LT.  1.0))  THEN 
fzdb  -  10.0  -  200.0  *  (zl  -  0.5)**4 
END  IF 

IF  ((Freq  .GT.  150.0) .AND. (Freq  .LE.  350.0) 

1  .AND. (zl  .GE.  1.0))  THEN 

fzdb  -  7.5  *  zl**( -13 . 3)  -  10.0 
END  IF 

IF  ((Freq  .GT.  350 .0) . AND. (zl  .GE.  1.0))  THEN 
fzdb  -  12.5  *  zl**( -8.0)  -  15.0 
END  IF 
ELSE 
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c  Calculate  evaporation  duct  height-gain  factor. 

z2  -  AMAXl(h  *  zfac,  hmin) 

IF  (Del  .GE.  10.25)  THEN 

c  Calculate  height-gain  for  del>-10.25  meters. 

IF  (z2  .GT.  zmax)  THEN 
fzdb  -  c5  *  (z2**c6)  +  c7 
ELSE 

fzdb  -  cl  *  AL0G(SIN(c2  *  (z2**c3)))  +  c4 
END  IF 
ELSE 

c  Calculate  height-gain  for  del<10.25  meters, 

fzdb  -  (cl  *  z2**c2)  +  (c3  *  z2**c4)  +  c5 
END  IF 
END  IF 

RETURN 

END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


Subroutine  insrt 


INSRT  inserts  (or  appends)  a  new  level  into  the  M-unit  profile.  It 
does  this  by  locating  the  new  height  relative  to  the  existing  pro¬ 
file  heights.  If  the  new  height  is  greater  than  the  top  level,  then 
append  a  new  level  for  the  new  height.  If  the  new  height  is  between 
two  levels,  then  Insert  a  new  level  for  the  new  height.  If  the  new 
height  is  equal  to  an  existing  level's  height,  do  not  add  a  new 
level  for  the  new  height. 


Variable : 
amu 
hmrs 


ipnt 


Description: 

Modified  refractivity  array,  M-units. 

Height  array,  meters,  each  element  corresponding  to 
the  like -number  amu  array  element. 

Number  of  levels  in  amu  and  hmrs. 

Height  of  new  level  to  be  added,  meters. 

Index  pointer  to  new  level. 


SUBROUTINE  insrt(amu,hrars , iq.hgt, ipnt) 

real*4  amu(32) ,hmrs(32) ,hgt 
integer*2  iq.ipnt 

DO  i-l.iq 
ilevel-i 

IF(ABS(hgt-hrars( ilevel) ) . LE . 0. 01)  go  to  1020 
IF(hmrs(Ilevel) .GT .hgt)  go  to  1030 
END  DO 


Hgt  >  amu(iq) 
iq-Iq+1 
ipnt-iq 

grdnt-0. 1181102 

amu(ipnt)-amu(iq-l)  +  (hgt-hmrs( iq- 1) )*grdnt 
hmr8( ipnt)-hgt 
go  to  1050 
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c 

c  Hgt  -  hmrs(ilevel) 

1020  continue 

ipnt-llevel 

c  amu(ipnt)-amu(llevel) 
hmrs(ipnt)-hgt 
go  to  1050 
c 

c  Hmrs(ilevel)  >  hgt  >  hmrs( 1 level -1) 

1030  continue 

c  Shift  levels  above  new  height  up  one 

DO  i-ilevel,iq 

1-iq  -  (i-ilevel) 
nmrs(j+l)-hmrs( j ) 
amu( j+l)-amu( j ) 

END  DO 

iq-iq+1 

ipnt-llevel 

grdnt-(amu(ipnt+l) -amu(ipnt-l) )/(hmrs(lpnt+l) -hmrs(ipnt-l)) 
amu(lpnt)-aniu(lpnt-l)  +  (hgt-hmrs(ipnt-l))*grdnt 
hmrs( lpnt)-hgt 
c  go  to  1050 
c 

1050  continue 
RETURN 
END 


c 

c  Subroutine  mprof 
c 

c  MPROF  modifies  the  M-unit  and  height  arrays  by  inserting  a  level  at 
c  the  antenna  height  using  straight  line  interpolation  (or  a  standard 
c  atmosphere  gradient)  to  calculate  its  M-unit  value.  The  new  profile 
c  is  then  used  to  locate  any  ducts  that  might  be  contained  in  the  pro- 
c  file.  If  the  bottom  of  tne  duct  is  below  the  EM  system  antenna 
c  height,  and  the  top  above  the  antenna  height,  then  a  critical  angle 
c  is  calculated  for  the  EM  system  in  the  surface -based  duct.  (It  is 
c  assumed  that  low-elevated  ducts  are  surface  ducts  if  the  EM  system  is 
c  in  the  duct.) 
c 

c  Variable: 
c  alphac 

c  amu 

c  antena 

c  antmu 

c  dcts 

c 
c 
c 

c  dMdh 

c  hbot 

c  htop 

c  height 

c  hmrs 

c 

c  lvlant 

c  lvltop 


Description: 

The  critical  penetration  angle  necessary  to  escape  duct 
An  array  of  M-unit  values 
EM  system  antenna  height 

M-unit  value  at  the  EM  system  antenna  height 
24  duct  parameter  array 

1, n  oottom  of  duct  'n' ,  meters 

2, n  top  of  duct  'n' ,  meters 

3, n  minimum  refractivity  of  duct  'n' ,  m-units 
M-unit  gradient  array 

Height  of  the  bottom  of  a  duct 
Height  of  the  top  of  a  duct 

Height  array  with  the  original  profile  heights 
Height  array  with  elements  corresponding  to  the  dMdh 
array  elements 
EM  system  antenna  level 

Maximum  number  of  layers  in  the  hmrs  array 
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Munits 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ndcts 
run  ax 
ntot 
rma 
sbdht 


M-unit  array  with  elements  corresponding  to  the  height 
array  elements 

The  number  of  ducts  stored  in  'dots' 

The  number  of  elements  in  the  height  and  Munit  arrays 
The  number  of  elements  in  the  dMdh  and  hmrs  arrays 
M-unit  value  at  the  minimum  on  the  duct  profile 
The  height  of  the  surface-based  duct 


Variables  not  listed  are  temporary  variables. 


SUBROUTINE  mprof (height , Munits , antena , nmax , ALPHAC , DMDH, HMRS , 
1  SBDHT, NTOT) 


c 

c 

real*4  alphac , amu( 32) , antena, dmdh( 32) ,hmrs(32) ,height(30) 
real*4  Mun its( 30 ), sbdht 
real*4  antmu.dcts ,hb ,ht , rma 
integer*2  lvlant , Ivltop , nmax, ntot 
integer*2  ndcts 
dimension  dcts(3,8) 
c 

Ivltop  -  nmax 
alphac  -0.0 
sbdht  -  0.0 

c 

c  Copy  height  and  m-unit  arrays, 

c 

Ivltop  -  nmax 
DO  i  -  1 ,  nmax 

hmrs(i)-height(i) 
amu(i)-Munits(i) 

END  DO 

Insert  new  level  at  the  antenna  height. 

call  insrt(amu,hmrs , Ivltop .antena, lvlant) 
antmu-amu( lvlant) 

Locate  all  major  ducts. 
ndcts-8 

call  due ts(amu, hmrs , Ivltop ,dcts .ndcts) 


c  Define  trapping  duct  parameters. 

IF(ndcts  .NE.  0)THEN 
DO  iduct-1 , ndcts 
hb-dcts( 1 , lduct) 
ht-dcts(2 , lduct) 
rma-dcts(3 , iduct) 

IF((antena  .GT.  hb)  .AND.  (antena  .LT.  ht))  go  to  1040 
IF(hb.lt.O.Ol)  go  to  1040 
END  DO 
END  IF 
c 

c  Antenna  not  inside  a  major  duct, 

go  to  1050 
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c  The  anCenna  is  inside  a  low-level  elevated  duct 

c  or  inside  a  surface -based  duct. 

1040  continue 
sbdht  -  ht 

alphac-1.0e-3*sqrt(2.0*(antmu-rma))  +  1.0e-5 
1050  continue 

Delete  all  levels  between  the  surface  and  the  antenna  level. 
DO  i  -  lvlant , lvltop 
j-i- (lvlant-2) 
hmrs(1 )-hmrs(i) 
a®u(j)-amu(i) 

END  DO 
lvltop-j 
lvlant-2 

c  Calculate  the  M-unlt  gradient  array. 

iend-lvltop-1 
DO  i  -  1,  lend 

delu-amu(i+l) -amu(i) 
delh-hmrs(i+l) -hmrs( i) 
dmdh(i)-l .0e-3*delu/delh 
END  DO 

dmdh<lvltop)-0 . 1181102e -  3 
c 

ntot  -  lvltop 

RETURN 

END 


c 

c 

c 

c 

c 

c 


Subroutine  opens t 

OPCNST  initializes  optical  region  constants. 

Variable:  Description: 

Effective  earth  radius,  (rk  *  6371),  km. 

Effective  earth  radius  *  1000. 

Aeth  *  2 

Dielectric  constant  of  sea-water,  epsilon. 

EM  system  frequency  in  MHz. 

Free-space  loss  constant  in  dB. 

RMS  wave  height  due  to  wind  in  m. 

Constant  for  subroutine  ruff, 

(hbar  *  2  *  PI  /  wavelength). 

EM  system  antenna  polarization: 

H  -  horizontal,  V  -  vertical,  C  -  circular 
Effective  earth  radius  factor. 

Real  part  of  the  square  of  the  index  of  refraction. 
Imaginary  part  of  the  square  of  the  index  of  refract. 
Constant  used  to  calculate  path- length  difference 
between  direct  and  sea-reflected  rays. 

Constant  (ae  *  2) . 

Wind  speed  in  kts . 


c 

ae 

c 

aeth 

c 

ae2 

c 

eps 

c 

freq 

c 

fsterm 

c 

hbar 

c 

hbfreq 

c 

c 

polar 

c 

c 

rk 

c 

rnreal 

c 

rnimag 

c 

thefac 

c 

c 

twoae 

c 

wind 

c 

c 
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SUBROUTINE  opcnst 


real*4  eps ,  sigma 

include  ' ffac . common' 
include  'envsys. common' 

fstenn  -  32.44  +  20.0  *  ALOGlO(freq) 
c  Exclusively  for  REF  subroutine 

IF  (polar  .NE.  "H")  THEN 

c  eps  is  the  permittivity  of  salt  water 

c  sigma  is  the  conductivity  of  salt  water 

IF  (freq  .LE.  1500.0)  THEN 
eps  -  80.0 
sigma  --4.3 

ELS El F  (freq  .LE.  3000.0)  THEN 

eps  -  80.0  -  0.00733  *  (freq  -  1500.0) 
sigma  -  4.3  +  0.00148  *  (freq  -  1500.0) 

ELSEIF  (freq  .LE.  10000.0)  THEN 

eps  -  69.0  -  0.00243  *  (freq  -  3000.0) 
sigma  -  6.52  +  0.001314  *  (freq  -  3000.0) 

ELSE 

eps  -  51.99 
sigma  -  15.718 
END  IF 

c  Define  the  real  and  imaginary  parts  of  the  square  of 

c  the  index  of  refraction  of  sea- water, 

rnreal  -  eps 

rnimag  -  (-18000.0)  *  sigma/freq 
END  IF 

c  Define  rms  wave -height  for  subroutine  RUFF 

hbar  -  0.0051  *  (0 . 5l4?7*wind)**2 
hbfreq  -  0.02094  *  freq  *  hbar 
c 

ae  -  rk  *  6371.0 
twoae  -  2.0  *  ae 
aeth  -  rk  *  6 . 371 
ae2  -  aeth  *  2.0 
thefac  -  freq  *  4.193E-5 
c 

RETURN 

END 
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c  Subroutine  opf£ac 
c 

c  OPFFAC  calculates  quantities  used  to  determine  the  pattern 
c  propagation  factor  (F)  in  the  optical  interference  region, 
c 

c  Variable: 
c  ae 

c  alpha 

c  angle 

c  beta 

c  divfac 

dr 


gamma 
psi 
rl 
r2 

range 
c  rmag 

c  ruf 

c  sinpsi 

c  twoae 

c 
c 

SUBROUTINE  opffac( gamma , range ,psi ,rl , r2 , rmag, ELANG , DPAT , DR) 

c 

real*4  angle,  beta,  divfac,  dpat,  dr,  elang,  gamma,  psi,  rl, 

1  r2,  range,  rmag,  ruf,  sinpsi 

c 

include  ' ffac . common' 
include  ' envsys . common' 

c 

patfac  -  1.0 

c  Calculate  direct  ray  launch  angle,  alpha, 

alpha  -  hdif/range  -  range/twoae 
angle  -  alpha 
elang  -  alpha 

c  Determine  antenna  pattern  factor  for  direct  ray  alpha, 

call  antpat (antype, alpha, antbwr.antelr .untfac, 

1  patrfac, angle, PATFAC) 

patd  -  patfac 
dpat  -  patfac 
beta  -  -  (gamma  +  psi) 
angle  -  beta 

c  Determine  antenna  pattern  factor  for  reflected  ray  beta, 

call  antpat (antype .alpha , antbwr , antelr , antfac , 

1  patrfac , angle , PATFAC) 

c  Determine  surface  roughness  coefficient, 

sinpsi  -  SIN(psi) 

call  ruff(hbar,  hbfreq ,  psi,  sinpsi,  RUF) 
c  Calculate  the  divergence  factor. 

divfac  -  1 . 0/(SQRT(l . 0  +  (2.0  *  rl  *  r2)/(ae  *  range  *  sinpsi))) 
dr  -  patfac  *  ruf  *  divfac  *  rmag 
c 

RETURN 

END 


Description: 

Effective  earth  radius,  km. 

Direct  ray  launch  angle,  radians. 

Angle  for  which  antenna  pattern  factor  desired. 
Reflected  ray  launch  angle,  radians. 

Divergence  factor. 

Constant  -  product  of  antenna  pattern  factor  for 
reflected  ray  *  divergence  factor  *  reflection 
coefficient  *  surface  roughness  factor. 

Earth's  interior  angle  (rl/ae). 

Grazing  angle,  radians. 

Reflection  point  range,  km. 

Reflection  point  range,  km. 

Total  ground  range  in  km. 

Magnitude  of  reflection  coefficient. 

Sea-surface  roughness  coefficient. 

Sin(psi) . 

2*ae 
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c 

c  Subroutine  oplimit 
c 

c  OPLIMIT  calculates  the  maximum  range  in  the  optical  region,  opmaxd, 
c  and  opmaxl  -  -20  LOG(F)  at  opmaxd,  where  F  is  the  pattern  propagation 
c  factor, 
c 

c  Variable: 

c  ae 

c  ae2 

c  alpha 

c  alphac 

c  al 

v  exloss 

c 

c  freq 

c  fsord 

c  gamma 

c  hdif 

c  horznl 

c  hi 

c  hip 

c  h2 

c  h2p 

c  opmaxd 

c  opmaxl 

c  pd 

c 

c  phi 

c  psilim 

c  psi 

c  r 

c  rl 

c  r2 

c  rk 

c  theta 

c 

c  thefac 

c  thnext 

c 
c 

SUBROUTINE  oplimit (OPMAXD, OPMAXL) 

c 

real*4  al,  dr,  fsqrd,  gamma,  halfpi,  hip,  h2p,  pd,  phi, 

1  pi,  psi,  psilim,  r,  rl,  r2 ,  rmag,  theta,  thnext 

c 

include  '  ffac . common' 

Include  ' envsys . common' 
c 

PI  -  3.14159 
halfpi  -  PI  /  2.0 
horznl  -  3.572  *  SQRT(rk  *  hi) 
psilim  -  0 . 01957/( f req*rk)**0 .33333 
c  If  both  ( erminals  are  in  the  duct  set  alphac  -  0.0 

IF  ((alphac  .GT.  0.0)  .AND.  (h2  ,LT.  sbdht))  alphac  -  0.0 
c  Initial  guess  for  rl  is  based  on  grazing  angle  limit  range, 

c  Use  ray  trace  equations  to  determine  rl  and  r2 . 

psi  -  psilim 

al  -  SQRT(psi**2  +  2 ,0e - 3*hl/ae) 
rl  -  (al  -  psi)*ae 
r2  -  rl 

IF  (h2  .GT.  hi) r2  -  r2  +  (SQRT(al**2  +  2 .0*ABS(hdif)/ae)  -  al)*ae 

r  -  rl  +  r2 

hip  -  hi  -  rl*rl/ae2 


Description: 

Effective  earth  radius,  km. 
ae  *  2000 

Direct  ray  launch  angle,  radians. 

Critical  angle  in  radians. 

An  angle  used  to  determine  rl(psllim). 

A  measure  of  how  much  of  the  antenna's  energy 
is  directed  toward  the  horizon,  dB. 

EM  system  frequency,  MHz. 

Square  of  the  pattern  propagation  factor,  F. 

Earth's  interior  angle. 

Difference  in  height  between  hi  and  h2. 

Tangent  ray  distance  for  height  hi,  km. 

Transmitter  height,  m. 

Effective  transmitter  height,  m. 

Receiver/target  height,  m. 

Effective  receiver/target  height,  m. 

Maximum  range  in  optical  region,  km. 

Propagation  factor  in  dB  at  opmaxd. 

Path-length  difference  between  direct  and 
sea- reflected  rays. 

Phase  lag  due  to  reflection  from  sea-surface. 

Grazing  angle  limit  to  optical  region. 

Grazing  angle  in  radians. 

Total  ground  range,  km. 

Reflection  point  range  (from  hi),  km. 

Reflection  point  range  (from  h2) ,  km. 

Effective  earth  radius  factor. 

Total  phase  difference  between  direct  and  sea-reflected 
rays  (pd)  and  phase-lag  due  to  reflection,  phi. 
Constant  used  to  calculate  path-lenght  difference. 

The  next  value  of  theta  to  be  determined. 
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h2p  -  h2  -  r2*r2/ae2 

call  ref(polar ,psi,RHAG,PHl) 

pd  -  thefac*hlp*h2p  /  r 

c 

c  Calculate  theta  based  on  grazing  angle  limit . 

c 

theta  -  phi  +  pd 

alpha  -  hdif/r  -  r/twoae 

IF  (alphac  .GT.  0.0)  THEN 

IF  ((alpha  .LT.  alphac)  .OR.  (pd  .GT.  halfpi))  THEN 
c 

c  Calculate  theta  based  on  range  obtained  from  alphac. 

c 

r  -  (SQRT(alphac**2  +  2 .0*ABS(hdif )/ae)  -  alphac)*ae 
call  opticf (polar , r , PD , PSI , THETA , FF) 

END  IF 
END  IF 

IF  ((alphac  .GT.  0.0)  .AND.  (pd  .GT.  halfpi);  THEN 
c 

c  If  theta>(2  Pi)  then  optical  limit  is  1st  peak 

c  with  theta  greater  than  theta(alphac) . 

IF  (theta  .GT.  6.28319)  THEN 

thnext  -  INT(theta/(2.0*PI)  +  1)*(2.0  *  PI) 
call  rliter(polar, thnext, R1,R2,R, PSI ,RMAG) 
theta  -  thnext 
?ND  IF 
ELSE 

Optical  limit  is  grazing  angle  limit  or  1/4  wavelength  limit. 
IF  ((pd  .GT.  halfpi)  .OR.  (psi  .NE.  pallia))  THEN 
Determine  theta  value  @  1/4  wavelength  limit,  (H  polar), 
thnext  -  1.5  *  PI 

call  rliterC'H", thnext ,R1 ,R2 ,R, PSI ,RMAG) 

IF  (polar  .NE,  "H")  THEN 

call  ref(polar ,psl,RMAC,PHI) 
theta  -  halfpi  +  phi 
ELSE 

theta  -  thnext 
END  IF 
END  IF 
END  IF 

IF  (ht  .GE.  hr)  THEN 
gamma  -  r2/ae 
ELSE 

gamma  -  rl/ae 
END  IF 

call  opffac( gamma , r ,psi , rl , r2 , rmag, ALPHA, PATD, DR) 
fsqrd  -  (patd*patd  +  dr*dr  +  ?. 0*dr*patd*COS( theta) ) 
c  Limit  fsqrd  to  prevent  runtime  errors  when  taking  LOG(fsqrd). 

IF  (fsqrd  ,LT.  1.0e-7)  fsqrd  -  1.0e-7 
cpmaxd  -  r 

opmaxl  -  *  10.0  *  ALOG10( fsqrd) 
exloss  -  -  20.0  *  ALOGlO(patd) 

RETURN 

END 
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c  Subroutine  opticf 
c 

c  Subroutine  OPTICF  calculates  the  total  phase  difference,  theta, 
c  between  direct  and  sea-reflected  ray  paths,  including  phase 
c  change  due  to  reflection  from  sea- surface.  It  then  uses  theta 
c  to  determine  the  value  of  the  pattern  propagation  factor,  F,  in 
c  the  optical  region,  and  returns  20Log(F). 
c 
c 

c  Variable: 
c  ae 

c  ae2 

c  aeth 

c  alpha 

c  dr 

c 
c 

c  epsr 

c  ff 

c  fprl 

c  frl 

c  fsqrd 

c  gamma 

c  hrp 

c  htp 

c  hi 

c  h2 

c  pd 

c 

c  phi 

c  psi 

c  r 

c  rl 

c  rlsqrd 

c  r2 

c  rr 

c  rmag 

c  t 

c  theta 

c 

c  thefac 

c  v 

c  w 

c 
c 

SUBROUTINE  opticf (plr , r , PD , PSI , THETA , FF) 
c 

real*4  dr,  epsr,  ff,  frl,  fprl,  fsqrd,  gamma,  hrp,  htp, 

1  phi,  psi,  r,  rl,  rlsqrd,  r2 ,  rmag,  rr,  t,  theta,  v,  w 

character*!  plr 
integer*2  jk 

c 

include  ' ffac .common' 
include  ' envsys . common' 
c 

rl  -  (hl/(hl  +  h2))*r 
t  -  1 .  b  *  r 

v  -  0.5  *  r  *  r  -  aeth  *  (hi  +  h2) 
w  -  aeth  *  r  *  hi 
epsr  -  0.050 
rr  **  2.0  *  epsr 
Jk  -  1 

C 


Description: 

Effective  earth  radius,  km. 

Ae*20U0,  km. 

Ae*1000,  km. 

Direct  ray  launch  angle,  radians. 

Product  of  divergence  factor,  surface  I'oughness 
coefficient,  reflection  coefficient  and  antenna 
pattern  factor  for  the  reflected  ray. 

Iteration  loop  range  tolerance,  km. 

Pattern  propagation  factor,  F,  in  dB. 

Value  of  the  derivative  of  che  cubic  equation  at  rl . 
Value  of  the  cubic  equation  for  a  given  rl. 

Square  of  the  pattern  propagation  factor. 

Earth's  interior  angle  (rl/ae)  in  radians. 

Effective  receiver/target  height,  m. 

Effective  transmitter  height,  m. 

The  transmitter  height,  di. 

The  receiver/target  height,  m. 

The  path-length  difference  between  direct  and  re¬ 
flected  rays  in  radians. 

Phase  lag  due  to  reflection  from  sea  surface,  rad. 
Grazing  angle  in  radians. 

Total  ground  range,  km. 

Reflection  point  range,  (from  xn.tr),  km. 

Square  of  trie  reflection  point  range. 

Reflection  point  range,  (from  revr/target) ,  km. 
Iteration  loop  variable  -  range  difference. 

Magnitude  of  reflection  coefficient. 

Iteration  loop  variable. 

Total  phase  lag  between  direct  and  sea-reflected 
rays,  in  radians.  (theta  -  pd  +  phi) 

Constant  used  to  calculate  theta. 

Iteration  loop  variable. 

Iteration  loop  variable. 
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c  WHILE  ((jk  .LT.  10) .AND. (abs(rr)  .GT.  epsr)) 
c 

c  Use  Newton-Raphson  iteration  to  solve  Kerr's  cubic  equation 

c  for  reflection  point  range  of  the  sea-reflected  ray.  (This 

c  equation  may  be  solved  explicitly  using  an  inverse  cosine.) 

c  The  Newton  iteration  works  best  if  hi  Is  less  than  h2. 

c 

DO  WHILE  ((jk  .LT.  10) . AND . (abs(rr )  .GT.  epsr)) 
jk  -  jk  +  1 
rlsqrd  -  rl  *  rl 

c  Kerr's  cubic  equation  for  reflection  point  range. 

frl  -  rl*rlsqrd  +  t*rlsqrd  +  v*rl  +  w 
c  Derivative  of  the  cubic  equation, 

fprl  -  3.0*rlsqrd  +  2.0*t*rl  +  v 
rr  -  frl/fprl 
rl  -  rl  -  rr 

IF  ((rl  .LT.  0.0) .OR. (rl  .GT.  r))  rl  -  r/2.0 
c  WEND 

END  DO 

r2  -  r  -  rl 

htp  -  hi  -  rl*rl/ae2 

hrp  -  h2  -  r2*r2/ae2 

psi  -  1.0e-3  *  htp  /  rl 

IF  (psi  .GT.  0.3)  psi  -  ATAN ( 1 . Oe - 3  *  htp/rl) 
call  ref(plr,psi,RMAG,PHI) 
pd  -  thefac*htp*hrp/r 
theta  -  pd  +  phi 
IF  (ht  . GE.  hr)  THEN 
gamma  -  r2/ae 
ELSE 

gamma  -  rl/ae 
END  IF 

call  opf f ac( gamma, r, psi , rl ,r2 , rmag, ALPHA, PATD, DR) 
fsqrd  -  patd*patd  +  dr*dr  +  2 . 0*dr*patd*COS ( theta) 
c  Limit  F-factor  to  -70  dB. 

IF  (fsqrd  .LT.  1.0e-7)  fsqrd  -  1.0e-7 
ff  -  -  10.0  *  ALOG10( fsqrd) 
c 

RETURN 

END 
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c 

c  Subroutine  push 
c 

PUSH  stores  elements  in  an  array  and  returns. 


Variable:  Description: 


array 

iq 

nq 

data 


iq  array  to  hold  data  elements 
Size  of  data  array 

Number  of  data  elements  stored  in  data  array 
The  data  element  to  be  stored 


SUBROUTINE  push(ARRAY, iq ,nq .data) 

real*4  data, array 
integer*2  iq.nq 
dimension  array(iq) 

Shift  array  elements  down  one 
do  1-iq, 2,-1 


array(i)-array(i- 1) 

END  DO 

Store  new  data  element  in  top  of  array 
array( l)-data 
nq-nq+1 

IF(nq  .GT.  iq)  nq  -  iq 

RETURN 

END 


Subroutine  rliter 


R1ITER  determines  a  reflection  point  range  'rl'  corresponding 
to  'rtheta'.  The  desired  reflection  point  range  is  determined  by 
a  Newton-Raphson  iteration  technique  to  vary  the  reflection  point 
point  range  until  the  correct  value  is  found. 


Vat iable : 
r 

rl 

r2 

f 

fl 

irlmda 

phi 

plr 

psi 

r 


Description : 

Distance,  or  range,  in  km. 

Distance  from  the  transmitting  antenna  to  reflection 
point  in  km. 

Distance  from  the  target/receiver  antenna  to  the 
reflection  point  in  km. 

Function  (Total  path  difference  between  direct  and 
sea-reflected  rays:  Theta)  used  in  iteration  loop. 

Finite  derivative  of  f. 

Iteration  loop  counter. 

Phase -lag  due  to  sea-surface  reflection  -  radians. 

EM  system  polarization  [H  -  horizontal,  V  —  vertical, 
C  -  circular] . 

Grazing  angle  in  radians . 

Range ,  in  Km , 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


rl  Distance  from  the  transmitting  antenna  to  reflection 

point  in  km. 

r2  Distance  from  the  target/receiver  antenna  to  the 

reflection  point  in  km. 

rmag  Magnitude  of  the  reflection  coefficient, 

rtheta  The  desired  value  of  theta. 


SUBROUTINE  rliter(plr , rthv  ca ,R1 ,R2 , R, PSI , RMAG) 


real*4  f,  f 1 ,  phi,  psi,  r,  rl ,  r2,  rmag,  rr,  rtheta 
character*!  plr 
integer*?,  irlmda 

include  ' ffac . common' 
include  'envsys .common' 

irlmda  -  0 
rr  -  rl 

WHILE  ( (abs(rr)  .GT.  0.001) .AND. (irlmda  .LT.  100)) 
(Equivalent  to:  100  IF  ((...). and. (...) )  THEN 

GOTO  100 

DO  WHILE  ((abs(rr)  .GT.  0.001) .AND. (irlmda  .LT.  100)) 

Calculate  phase  difference,  theta,  corresponding  to 
reflection  point  range  rl.  Then  use  finite  derivative 
method  to  iterate  to  the  range  where  theta  is  equal  to 
the  target  value:  rtheta. 

call  gtheta(plr,rl,R,F,R2, PSI, RMAG) 

call  gtheta(plr , rl+0 . 001 , R , FI , R2 , PSI , RMAG) 

fp  -  (fl  -  f)  /  0.001 

rr  -  (rtheta  -  f)  /  fp 

irlmda  -  irlmda  +  1 

IF  (rr  .GT.  -rl)  THEN 

IF  (rr  +  rl  .LE.  horznl)  THEN 
rl  -  rl  +  rr 

ELSE 

rl  -  (rl+horznl)/2 .0 
END  IF 
ELSE 

rl  -  r 1/2.0 
END  IF 
END  DO 
WEND 
RETURN 
END 
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c 

c  Subroutine  ref 
c 

c  Subroutine  REF  returns  the  magnitude  and  phase  lag  of  the  reflection 
c  coefficient  for  reflection  from  the  (smooth)  sea  surface.  These 
c  quantities  are  calculated  as  a  function  of  the  grazing  angle  psi. 
c  The  complex  square  roots  are  done  by  separating  the  complex  variables 
c  into  their  real  and  imaginary  parts.  Ho  complex  function  calls  are 
c  used, 
c 
c 

c  VARIABLE: 
c  rnreal 

c 

c  rnimag 

c 
c 

c  phi 

c  plr 

c 

c  psi 

c  mag 

sinpsi 


DESCRIPTION: 

Real  part  of  the  square  of  the  index  of  refraction, 
(the  dielectric  constant  of  sea-water). 

Imaginary  part  of  the  square  of  the  index  of 
refraction  (the  conductivity  of  sea  water 
times  the  wavelength  times  other  constants). 

Phase  change  (lag)  in  radians. 

EM  system  antenna  polarization:  H  -  horizontal; 

V  -  vertical;  C  -  Circular. 

Crazing  angle  in  radians. 

Magnitude  of  the  reflection  coefficient. 

SIN(psi) . 


c 

c 


c 


c 


c 


c 


c 


various  All  variables  not  listed  above  are  temporary. 
SUBROUTINE  ref (plr,  psi,  RMAC ,  PHI) 

real*4  angrt,  at,  bt,  ct,  dt,  phi,  phiv,  pi,  psi, 

1  rev,  rmag,  rmagrt,  rtimag,  rtreal, 

2  rvimag,  rvreal,  rx,  sinpsi,  x,  y 
character*l  plr 

Include  ' ffac . common' 
include  'envsys .common' 


PI  -  3.14159 

Define  RMAG,  PHI  for  horizontal  polarization, 
rmag  -  1.0 
phi  -  PI 

IF  (plr  .NE.  "H")  THEN 

Calculate  RMAG,  PHI  for  vertical  polarization, 
sinpsi  -  SIN(psi) 

Y  -  rnimag 

X  -  rnreal  -  C0S(psi)**2 

rmagrt  -  (x*x  +  y*y)  **  0,25 

angrt  -  ATAN(y/x)  /  2.0 

rtreal  -  rmagrt  *  COS(angrt) 

rtimag  -  rmagrt  *  SIN(angrt) 

at  »  rnreal  *  sinpsi  -  rtreal 

ct  -  rnreal  *  sinpsi  +  rtreal 

bt  -  rnimag  *  sinpsi  -  rtimag 

dt  -  rnimag  *  sinpsi  +  rtimag 

rvreal  -  (at*ct  +  bt*dt)  /  (ct**2  +  dt**2) 

rvimag  -  (bt*ct  -  at*dt)  /  (ct**2  +  dt**2) 

rev  -  SQRT(rvreal+*2  +  rvimag**2) 

IF  (rvreal  .NE.  0.0)  THEN 

?hiv  -  ATAN( rvimag/ rvreal) 

F  (rvreal  ,LT,  0.0)  phiv  -  phiv  +  PI 
ELSE 


IF  (rvimag  .LT.  0.0) 
IF  (rvimag  CT,  0.0) 
IF  (rvimag  .EQ.  0.0) 
END  IF 


phiv  -  -PI  /  2.0 
phiv  -  PI  /  2,0 
phiv  -  0,0 


A-30 


phiv  -  -phlv 

IF  (phlv  .LT.  0.0)  phlv  -  phiv  •*-  2.0*PI 
rmag  -  rev 
phi  -  phiv 

IF  (plr  .EQ.  "C")THEN 

c  Calculate  RMAG ,  PHI  for  circular  polarization. 

rx  -  SQRT (1.0  +  rcv**2  +  2.0*rcv  *  C0S(PI  -  phiv)) 
rmag  -  rx/2 . 0 

a  -  rev  *  SIN(phiv  +  PI)  /  rx 
a  -  ATAN(  a/SQRT ( 1  -  a*a)  ) 
phi  -  PI  -  a 
phi  -  *phl 

IF  (phi  .LT.  0.0)  phi  -  phi  +  2.0*PI 
END  IF 
END  IF 
RETURN 
END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Subroutine  ruff 

Subroutine  RUFF  returns  the  sea-surface  roughness  correction  for 
the  magnitude  of  the  sea-reflected  ray. 


VARIABLE: 

hbar 

hbfreq 

hfpai 

psl 

slnpsi 

rufeo 


DESCRIPTION: 

rms  wave  height  in  meters. 

( 2*PI*hbar )/wavelength . 

(hbar*psi) /wavelength. 

Crazing  angle  in  radians. 

SIN(psi). 

Sea-surface  roughness  coefficient. 


SUBROUTINE  ruff (hbar,  hbfreq,  psi,  sinpsi,  RUFC0) 
real*4  hbar,  hbfreq,  hfpsi,  psi,  rufeo,  slnpsi 


rufeo  -  1.0 

IF  (hbar  .NE.  0,0)  THEN 

hfpsi  -  hbfreq  *  psl  *  0.159155 
IF  (hfpsi  .LE.  0.11)  THEN 

rufeo  -  EXP((-2,0)  *  (hbfreq*sinpsi)**2) 

ELSEIF  (hfpsi  .LE,  0.26)  THEN 

rufeo  -  0.5018913  -  SQRT(0 . 2090248  -  (hfpsi  -  0.55189)**2) 
ELSE 

rufeo  -0.15 
END  IF 
END  IF 
RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Subroutine  sysfil 

SYSFIL  list  available  system  files  and  allows  the  user  to  select  an 
EM  system  file. 


Variable : 
antena 

antype 


bwidth 

elevat 

filename 

freq 

polar 


rectar 


Description: 

Height  of  EM  system  antenna  in  m. 
Antenna  type: 

0  -  omnidirectional 
S  -  sin(x)/x 
C  -  cosecant -squared 
H  -  generic  height-finder 
Antenna  beam  width  in  degrees. 
Antenna  elevation  angle  in  degrees. 
Name  of  System  file. 

EM  system  frequency  in  KHz. 

Antenna  polarization: 

H  -  horizontal 
V  -  vertical 
C  -  circular 

Receiver/target  height,  m. 


SUBROUTINE  sysfil (freq,  antena,  rectar,  polar,  antype,  bwidth, 
1  elevat) 


c 


c 


c 


real*4  antena,  bwidth,  elevat,  freq 
integer*2  ZR,  ZW 
character*l  antype,  polar 
character*20  filename 

Initialize  read  and  write  channels. 
ZR  -  5 
ZW  -  6 


call  system  ('Is  (SJ*  1>&2 ' //char (0) ) 
write  (ZW, '(//, "Enter  input  file  name:  ",$)') 
read  (ZR,'(al2)')  filename 
open  (lO.FILE-filename) 


read  (10,  '(flO.l)’)  freq 
read  (10,  '(flO.l)')  antena 
read  (10,  '(flO.l)')  rectar 
read  (10,  '(al)')  polar 
read  (10,  '(al)')  antype 
read  (10,  '(flO.l)')  bwidth 
read  (10,  '(flO.l)')  elevat 
close(10) 


!  radar  freq 
!  radar  antenna  ht 
!  receiver/target  height 
!  antenna  polarization 
!  antenna  type 
!  vert,  beam  width 
!  ant.  elev.  angle 


RETURN 

END 
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c 

c  Subroutine  sysinp 
c 

c  Subroutine  SYSINP  prompts  the  user  for  EM  system  parameters  and  re- 
c  turns.  System  parameters  can  be  entered  from  the  keyboard  or  from  a 
c  file.  If  the  system  Is  entered  from  the  keyboard  it  can  be  stored  in 
c  a  file, 
c 

c  Variable : 
c  a type 

c 
c 
c 
c 

c  beam 

c  elang 

c  fmhz 

c  filename 
c  rectar 

c  plr 

c 
c 
c 

c  rectar 

c  xmtr 

c 
c 

SUBROUTINE  sysinp (FMHZ,  XMTR,  RECTAR,  PLR.  ATYPE,  BEAM,  ELANG) 
c 
c 

real*4  beam,  elang,  fmhz,  rectar,  xmtr 
character*l  atype,  dummy,  plr 
character*20  filename 
integer*2  ZW,  ZR 

c 

c  Specify  the  read  (5)  and  write  (6)  channel  numbers. 

ZW  -  6 
ZR  -  5 
c 

c  Enter  EM  system  parameters  from  file  or  keyboard. 

write (ZW, ' ("You  may  enter  EM  system  from  a  file  or  keyboard.")') 
wrlte(ZW ,'( "Enter  EM  system  data  from  file?  (yes  or  no)  ",$)') 
read(zr , ' (Al) ' )dumray 

IF  ( ( dummy (1:1)  .eq.  'y' ) .or. ( dummy ( 1 : 1 )  ,eq.  'Y'))  THEN 
c  Enter  EM  system  data  from  a  file. 

call  sysfil (fmhz , xmtr , rectar , plr , atype .beam , elang) 

ELSE 

c  Enter  the  EM  system  parameters  from  keyboard. 

wrlte(ZW, ' ("Enter  EM  System  Parameters:  ")') 
c 

c  Initialize  EM  system  variables, 

fmhz  -  5600.0 
xmtr  -  25.0 
rectar  -  25.0 
plr  -  "H" 
atype  -  "0" 
beam  -  0.0 
elang  -  0.0 
c 

vrite(ZW, 1010) 

1010  format( ' Enter  frequency  in  MHz  (100  to  20,000)  ' ,$) 
read(ZR,*)  fmhz 

IF  (fmhz  ,LT.  100.0)  fmhz  -  100.0 
IF  (fmhz  .GT.  20000.0)  fmhz  -  20000,0 


Description: 

Antenna  type : 

0  ~  omnidirectional 
S  -  sin(x)/x 
C  -  cosecant- squared 
H  -  height- finder 
Beam  width  in  degrees . 

Antenna  pointing  (elevation)  angle  in  degrees. 
EM  system  frequency  in  MHz. 

EM  system  filename. 

Receiver/radar  target  height,  m. 

Antenna  polarization: 

H  -  horizontal 
V  -  vertical 
C  -  circular 

Receiver/target  height,  m. 

EM  system  antenna  height,  m. 
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write(ZW, 1015) 

1015  format( ' Enter  transmitter  height  In  meters  (1  to  100)  ',$) 
read(ZR,*)  xmtr 
IF  (xmtr  .LT.  1.0)  xmtr  -  1.0 
IF  (xmtr  .GT.  100.0)  xmtr  -  100.0 
c 

write (ZW, 1020) 

1020  format( ' Enter  receiver/target  height  in  meters  (1  to  30000)  ',$) 
read(ZR,*)  rectar 
IF  (rectar  .LT.  1,0)  rectar  -  1.0 
IF  (rectar  .GT.  30000.0)  rectar  -  30000.0 
c 

wrlte(ZW, 1025) 

1025  format( ' Enter  EM  system  polarization  (H,  V,  C)  ',$) 
read(ZR, ' (Al) ' )  plr 

IF  ((plr  .EQ.  "c")  .OR.  (plr  .EQ.  MCB))  plr  -  "C" 

IF  ((plr  .EQ.  "v")  .OR.  (plr  .EQ.  "V"))  plr  -  "V" 

IF  ((plr  .NE.  "V")  .AND.  (plr  .NE.  "C"))  plr  -  "H" 
c 

wrlte(ZW, 1030) 

1030  format( ' Enter  antenna  type  -  options  are:  Omnidirectional,  ' 

1  /,'Sln(x)/x,  Cosecant- squared.  Height-finder  (O,  S,  C,  H))  ' ,$) 

read(ZR, ' (Al) ' )  dummy 

IF  ((dummy  .EQ.  "o")  .OR.  (dummy  .EQ.  "0"))  atype  -  B0" 

IF  ((dummy  .EQ.  "s")  .OR.  (dummy  .EQ.  "S"))  atype  -  "S" 

IF  ((dummy  .EQ.  "c")  .OR.  (dummy  .EQ.  "C"))  atype  -  "C" 

IF  ((dummy  .EQ.  "h")  .OR.  (dummy  .EQ.  "H"))  atype  -  "H" 

IF  ((atype  .NE.  "S")  .AND.  (atype  .NE.  "H")  .AND. 

1  (atype  .NE.  "C"))  atype  -  "0" 

c 

beam  -  0.0 
elang  -  0.0 

IF(atype  .NE.  "O")  THEN 
wrlte(ZW, 1035) 

1035  format (' Enter  antenna  beam  width  In  degrees  (>0.0  to  45)  ',$) 
read(ZR ,  *)  beam 
IF  (beam  ,LE.  0.0)  beam  -  0.10 
IF  (beam  .GT.  45.0)  beam  -  45.0 
c 

write(ZW, 1040) 

1040  format( ' Enter  antenna  elevation  angle  in  degrees  (-10.0  to  10.0)' 
1  /,'(0  Is  normal)  ',$) 

read(ZR,*)  elang 

IF  (elang  ,LT.  -10.0)  elang  -  -10.0 
IF  (elang  .GT.  10.0)  elang  -  10.0 
END  IF 

wrlte(ZW, ' ("Do  you  wish  to  store  this  EM  system  In  a  file?", 

1  *’  (yes  or  no)  ",$)') 

read(zr , ' (A) ' )dummy 

IF  ((dummy(l:l)  .eq.  ' y' ) .or . (dummy (1 : 1)  .eq.  'Y'))  THEN 
write  (ZW,'("  Current  System  Flies:  ")') 
call  system  ('Is  ( S ) *  1>&2 ' //char(0) ) 
write  (ZW,  1045) 

1045  format( "Enter  file  name  (First  letter  MUST  be  S):  ",$) 
read  (ZR, ' (al2) ' )  filename 
open  (10 , FILE-filename) 

c 

c  Write  frequency  and  antenna  heights  In  file, 

write (10,  '(flO.l)')  frahz 
write(10,  '(flO.l)')  xmtr 
wrlte(10,  '(flO.l)')  rectar 
c  Write  the  antenna  characteristics  in  file. 

write(10,  '(al)')  plr 
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write(10,  '(al)')  atype 
write(10,  '(flO.l)')  bean 
write(10,  '(flO.l)')  elang 
c  close  file 

close(lO) 

END  IF 
END  IF 
c 

RETURN 

END 


c 

c  Subroutine  tropo 
c 

c  Subroutine  TROPO  returns  the  troposcatter  loss  for  a  given  range, 
c  Troposcatter  loss  is  based  on  models  by  Yeh  with  a  frequency- gain 
c  correction  term,  HO,  from  National  Bureau  of  Standards  Document 
c  NBS  101.  Frequency  gain  factor  gives  additional  loss  for  low 
c  frequency,  low-sited  antennas, 
c 


c  VARIABLE: 

c  ae 

c  exloss 

c  fcubed 

c  horizn 

c  hi 

c  h2 

c  hO 

c  r 

c  rnsterm 

c  rnsubs 

c  rone 

c  rtwo 

c  tfac 

c  tloss 

c  tsubO 

c  tsubl 

c  tsub2 

c  ttot 

c 

c  various 

c 

c 


DESCRIPTION: 

Effective  earth  radius  in  kilometers 

Antenna  gain  for  lowest  optical  region  ray  in  dB 

Frequency  cubed 

Horizon  range  in  kilometers 

Transmitter  height  in  meters 

Radar  target/receiver  height  in  meters 

Frequency  gain  factor  in  dB 

Ground  range  in  km 

Constant  involving  Surface  Modified  refractivity 

Modified  refractivity  value  at  the  sea  surface 

4*PI*hl*ttot/wavelength 

4*PI*h2*ttot/wavelength 

Troposcatter  region  constant 

Troposcatter  loss  in  dB 

Angle,  theta  sub  0,  associated  with  total  range  r 

Angle,  theta  sub  1,  associated  with  horizon  range  rl 

Angle,  theta  sub  2,  associated  with  horizon  range  r2 

Scattering  angle,  (theta  in  NBS  101) 

All  variables  not  listed  above  are  temporary.  Most 
variables  use  the  names  given  in  NBS  101. 


SUBROUTINE  tropo(r , tloss) 


real*4  chi,  csubl,  csub2 ,  delhO,  etas,  horizn,  hsubO, 

1  hO,  hOrl,  h0r2,  q,  r,  rnsterm,  rnsubs,  rone,  rtwo, 

2  s,  tfac,  tloss,  tsubO,  tsubl,  tsub2,  ttot,  zeta 

include  ' ffac . common' 
include  ' envsys .common' 

rnsubs  —  Munits(l) 
tfac  -  0.08984  /  rk 

horizn  -  3.572  *  (  SQRT ( rk  *  hi)  +  SQRT ( rk  *  h2)  ) 
tsubO  -  r  /  ae 
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tsubl  -  SQRT(hl  *  ae/500.0)  /  ae 
tsub2  -  SQRT(h2  *  ae/500.0)  /ae 
ttot  -  CsubO  -  tsubl  -  Csub2 

zeta  -  ttot/2.0  +  tsubl  +  (hi  -  h2)  /  (1000. 0*r) 

chi  -  ttot/2.0  +  tsub2  +  (h2  -  hi)  /  (1000. 0*r) 

rone  -  hi  *  0.0419  *  freq  *  ttot 

rtwo  -  h2  *  0.0419  *  freq  *  ttot 

IF  (rone  ,LT.  0.1)  rone  -  0.1 

IF  (rtwo  .LT.  0.1)  rtwo  -  0.1 

s  -  zeta  /  chi 

IF  (s  .GT.  10.0)  s  -  10.0 

IF  (s  .LT.  0.1)  s  -  0.1 

q  -  rtwo  /  (s  *  rone) 

IF  (q  .GT.  10.0)  q  -  10.0 

IF  (g  .LT.  0.1)  q  -  0.1 

hsubO  -  s  *  r  *  ttot  /  (1.0  +  s)**2 

etas  -  0.5696*hsub0  *  (1.0  +  rnsterra*EXP(-3.8e-6  *  hsub0**6)) 

IF  (etas  .GT.  5.0)  etas  -  5.0 

IF  (etas  .LT.  0.01)  etas  -  0.01 

csubl  -  16.3  +  13.3*etas 

csub2  -  0.4  +  0.16*etas 

hOrl  -  csubl  *  (rone  +  csub2)**( -1 . 333) 

h0r2  -  csubl  *  (rtwo  +  csub2)**( -1 . 333) 

hO  -  (hOrl  +  h0r2)  /  2.0 

delhO  -  1.13  *  (0.6  -  ALOGlO(etas) )  *  ALOG(s)  *  ALOG(q) 

IF  (delhO  .GT.  hO)  THEN 
hO  -  2 . 0*h0 
ELSE 

hO  -  hO  +  delhO 
END  IF 

IF  (hO  .LT.  0.0)  hO  -  0.0 

tloss  -  114.9  +  tfac*(r-horizn)  +  10.0*ALOG10(r*r*freq**3) 
tloss  «  tloss  -  rnsubs+0.2  +  hO  +  exloss 

RETURN 

END 
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